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I .  INTRODUCTION 


In  this  grant  study  of  heat  conduction  in  composite  materials, 
several  phases  of  investigation  were  started  in  different  stages  of 
the  grant  period.  These  phases  are  individually  self-contained  from  a 
technical  viewpoint,  but  are  related  to  one  another  in  an  overall  sense. 
They  are  either  sequential  or  proceed  in  parallel  paths  leading  to  a 
common  objective. 

For  steady-state  applications  and  for  the  purpose  of  acquiring  a 
basic  data  bank,  a  methodology  to  analyze  the  effective  transverse 
thermal  conductivities  of  composites  consisting  of  uni-directional 
fibers  was  developed  in  Phase  I  (AFWAL-TR-80-3012,  Han  and  Cosner). 
Extensive  calculations  were  performed  to  obtain  accurate  values  for 
the  transverse  conductivities  of  composites  with  various  packing 
densities,  dispersion  patterns  and  tow-matrix  conductivity  ratios. 

The  large  amount  of  data-heretofore  unavailable,  but  generated  in 
Phase  I  -  was  used  subsequently  by  Zimmerman  in  a  Phase  II  report 
(AFWAL-TR-80-3155,  R.  H.  Zimmerman)  in  which  the  then-existing  predictive 
schemes  were  examined  in  detail.  Using  the  accurate  data  from  Phase  I 
as  a  base,  prior  simplified  predictive  equations  were  modified  to 
extend  their  ranges  of  validity.  Above  all,  a  unified  approach  to 
predict  the  effective  transverse  (to  fibers)  conductivities  was 
achieved  with  accuracy  of  5  to  10  per  cent.  Zimmerman's  work  was 
significant  in  Its  comprehensiveness  and  definitiveness  of  its 
conclusions. 
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As  thermal  transients  are  a  crucial  concern  in  the  analysis  of 
thermal  stresses  and  strains  -  such  as  delamination  -  investigations 
on  thermal  transients  were  initiated  in  this  grant  study,  and  they 
were  pursued  in  several  tasks.  One  task  is  concerned  with  the 
experimental  determination  of  directional  thermal  conductivities 
and  thermal  diffusivities  of  fibrous  materials.  The  purposes  of  this 
task  are  two-fold:  (1)  to  evolve  a  simple  and  reliable  experimental 
technique,  and  (2)  to  acquire  basic  data  for  graphite/epoxy  composites, 
with  future  extensions  to  other  (metallic)  composites.  This  phase  of 
work  will  be  documented  in  a  separate  report  shortly. 

A  second  task  is  concerned  with  the  development  of  a  simplified 
method  of  calculating  transient  temperature  distributions  in  composites 
with  problems  involving  heat  flows  along  the  uni -directional  fiber  as 
a  starting  point.  This  task  is  motivated  by  the  fact  that  rigorous 
solutions  of  the  governing  diffusion  equations  for  different  media 

in  composites  are  conceptually  speaking  possible,  but  practically 
speaking  prohibitive.  The  methodology  developed  is  that  of  a  heat- 
balance  Integral  method  -  a  concept  whose  origin  lies  in  the  work  of 
von  Karman  on  the  analysis  of  boundary  layer  flows. 

Transient  heat  flow  analyses  in  composite  bodies  are  not  new;  they 
date  back  to  one  hundred  years.  The  early  work  and  the  majority 
of  current  efforts  are  largely  patterned  after  classical  methods  of 
analysis  for  simplified  geometries.  Even  for  simplified  geometries, 
the  mathematical  complexities  are  quite  forbidding.  For  example,  the 
use  of  the  Laplace  transform  method  to  two-phase  problems  often  results 
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In  solving  a  non-linear  algebraic  equation  for  the  transform  parameter. 
Numerical  Implementation  of  Its  solution  necessitates  various  kinds  of 
perturbation  analyses,  and  the  amount  of  labor  required  renders  these 
conceptually  elegant  and  exact  method  much  less  elegant  and  less  exact, 
when  translated  Into  a  practical  level.  Hence,  a  need  exists  for  a 
more  pragmatic  approach  which  Is  more  accessible  to  design  engineers 
confronted  with  thermal  problems  of  composites.  The  contents  of  this 
report  are  a  first  step  to  achieve  the  goal. 

Presented  In  this  phase  III  report  is  the  development  of  a  heat- 
balance  integral  method  for  the  analysis  of  thermal  transients  in 
composites  with  two  different  media.  As  an  initial  attempt,  the 
integral -method  is  applied  to  heat  flow  problems  in  composites 
comprised  of  two  different  laminated  materials,  and  in  composites  with 
uni-directional  fibers  dispersed  evenly  in  a  binding  matrix.  Heat 
flow  takes  place  along  the  fiber  axial  direction  or  in  the  plane  of 
laminations  due  to  a  step-rise  in  the  surface  temperature  or  a  step- 
rise  of  the  surface  heat  flux.  The  latter  is  the  most  commonly 
encountered  boundary  condition  and  is  the  most  practical  one  to 
consider.  Problems  with  more  involved  boundary  conditions  can  be 
constructed  on  the  basis  of  a  combination  of  these  two  fundamental 
solutions. 

It  is,  of  course,  recognized  that  the  heat-balance  integral  method 
is  intrinsically  approximate  but  with  increasing  accuracy  as  the 
trial  functions  -  to  be  explained  later  in  this  report  -  satisfy  more 
secondary  boundary  conditions.  Therefore,  in  order  to  establish  and 
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confirm  Its  usefulness,  it  is  necessary  to  compare  with  available  exact 
solutions  for  the  problems  considered.  Consequently,  exact  solutions 
of  tiie  types  of  problems  described  earlier  form  the  core  objective  of 
another  task  of  the  grant  study.  The  rigorous  solution  is  either  in 
the  form  of  analytical  developments  or  in  the  form  of  numerical 
solutions  with  both  approaches  revealing  the  micro-structural 
temperature  distributions.  It  Is  to  be  emphasized  that  the  detailed 
computations  of  micro-structure  temperature-time  histories  are  not  an 
end  by  Itself.  They  serve  a  purpose  of  identifying  where  the 
heat-balance  integral -method  may  be  deficient  and  therefore  requires 
further  refinements.  A  report  outlining  this  particular  phase  of  work 
should  be  completed  soon  after  this  report  is  available. 
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II.  TWO  QUARTER- INFINITE  REGIONS 


As  an  introductory  Illustration  of  approximate  analyses  of  heat 
flow  in  composites,  consider  the  problem  of  two  quarter- infinite  regions 
of  two  different  materials,  1  and  2,  occupying  the  right-hand  side  of  the 
y-axis  in  a  x-y  space  (See  Figure  II-l).  In  the  first  quadrant, 
material  1  has  physical  properties  designated  by  p^,  and  Cj  and  in 
the  fourth  quadrant  material  2  has  properties  of  k^,  p2,  and  c2< 

With  both  regions  initially  at  a  uniform  temperature  (considered 
zero),  it  is  desired  to  determine  the  temperature  responses  in  these 
two  regions  when  the  temperature  on  the  surface  x  =  0  is  impulsively 
raised  to  a  finite  value  or  is  given  an  impulsive  heat  flux.  The  key 
feature  of  the  problem  is  of  course  the  interface  at  y  =  0  across 
which  inter-region  heat  conduction  affects  the  respective  temperature 
distributions  in  the  vicinity  of  the  interface.  Although  this  particular 
heat  conduction  problem  has  been  elegantly  studied  by  Grimado  (Quarterly 
of  Applied  Mathematics  p.379,  Jan.  1954),  the  results  obtained  by  him 
are  not  readily  usable  from  an  engineering  viewpoint.  With  a  minimum 
of  complications,  this  example  problem  serves  to  illustrate  the  essential 
features  of  the  heat-balance  integral  method,  which  are  applied  to 
problems  with  more  practical  interest  in  Sections  III  and  IV. 
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II.l  CONSTANT  SURFACE-TEMPERATURE  SOLUTIONS 


In  each  region  away  from  the  interface,  the  temperature  response 
is  asymptotically  identical  to  that  in  a  semi-infinite  region.  These 
asymptotic  solutions,  valid  for  |y|  -*■  «,  are: 


(II-l) 


(11*2) 


It  is  of  interest  to  note  that  it  is  possible  to  have  two  materials  with 
equal  diffusivity  (o]  a  ct2)  such  that  Equations  II-l  and  II-2  are 
identical.  In  this  case,  no  interface  distortion  would  result.  Our 
interest  is  of  course  the  general  case  when  f  a2- 

When  ct1  t  a2,  the  temperature  responses  given  by  Equations  II-l 
and  I 1-2  cannot  be  valid  in  a  region  near  the  interface  y  =  0.  The 
temperature  distribution  is  in  a  transition  stage  from  T^  at  y  =  « 
through  the  interface  and  then  to  T2<B  at  y  »  -®.  Figure  II. 2a 
Illustrates  such  a  distribution.  The  purpose  of  an  approximate  analysis 
is  to  determine  the  transitional  temperature  distribution. 
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Shown  In  Figure  II. 2b  Is  a  cross-sectional  view  of  the 
temperature  surface  In  Figure  II .2a.  The  temperature  distribution 
In  each  region  Is  conceptualized  Into  two  zones:  (1)  a  boundary 
zone  of  width  A  and  (ii)  the  rest  of  the  region.  In  the  boundary  zone, 
temperature  varies  from  an  interface  temperature  to  the  fully 
established  value  corresponding  to  either  one  of  Equations  II-l  and 
II-2.  In  other  words,  the  distribution  Is  simplified  by  the 
discretized  quantity  A^  or  A2  such  that  only  near  the  interface 
does  the  temperature  undergo  a  transition  from  one  established 
value  to  another. 

Suggested  in  Figure  II. 2b  are  therefore  two  boundary  zones 
A-|  and  Ag.  and  within  them,  the  temperature  distributions  for 
materials  1  and  2  are  given  by  the  following: 


W  (T,_-  T.lFly,/^) 

( 1 1-3 ) 

T2  *  T1  *  <Tj_  -  T,)F(y 2/42) 

( I 1-4) 

Note  that  the  same  functional  combination  is  used  on  either  side 
of  the  interface.  Equations  I 1-3  and  I 1-4  are  of  course  limited 
to  the  boundary  zones;  outside,  the  temperatures  ^  and  T2  are 
equal  to  and  T ^  directly.  The  latter  requirement  leads  to 
the  boundary  condition  on  F  at  y1  *  A1  or  y2  3  A2: 

F(1 )  -  1  ( I 1-5) 
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In  Equations  I 1-3  and  II-4,  the  common  interface  temperature 
is  denoted  by  y  which  both  and  Tg  must  assume  when  y^  *  y2  =  0. 

This  then  requires: 

F(0)  =  0  ( I 1-6) 

Finally  of  course  is  the  equality  of  the  heat  flux  at  the 
interface,  expressed  by: 

^Oyay,)  +  k2(3T2/3y2)  =  0 

Differentiating  Equations  II-3and  II-4  and  putting  y1  *  y2  =  0, 
there  results: 

T1  -  t(k1M1)T1.c8+  (k2M2)T2.ool/[(k1M1)  +  (fcjMg)]  (II-7) 


Iljy,! 1, I** 


1 1-2  THE  HEAT-BALANCE  INTEGRAL-METHOD 


The  governing  differential  equations  for  both  regions  are  of 
course  the  diffusion  equations: 


(3^/38)  = 


(3T2/36)  * 


(11-8) 

( 1 1-9) 


Since  the  above  equations  are  satisfied  by  T^^  and  T2-oo  defined  by 
Equations  II-l  and  I 1-2  for  these  two  separate  regions  (each  in  a 
semi-infinite  domain).  Equations  1 1-8  and  1 1-9  can  also  be  written 
in  terms  of  the  respective  differences  (T^  -  Tl  oo)  and  (T2  -  T2  oo), 
thus  resulting  in: 


-  T,j  •  a/tT,..  -  T,] 


(11-10) 


-  V  *  «/lT2~  ■  V 


(11-11) 


As  the  central  idea  of  the  heat-balance  integral -method  is 
that  a  differential  equation  Is  satisfied  on  the  average  while  the 
boundary  conditions  are  fulfilled  through  appropriate  functional 
forms.  Hence,  by  Integrating  Equation  11-10  from  y1  *  0  to 
y^  ■  the  result  is: 


12 


Vd*i  - 


1  -°° 


(11-12) 


A  similar  equation  is  obtained  starting  from  Equation  11-11 ;  the 
resulting  equation  is  obtainable  from  Equation  11-12  by  changing 
the  subscript  1  to  2. 

Now,  the  distribution  function  F  can  be  chosen  to  conform 
with  the  requirements  enumerated  before.  A  popular  choice  in 
viscous  boundary  layer  analyses  is  the  following: 


3/y.\  i  /y,y 
F(W  =  2\A^)  "  2\A^) 


(11-13) 


In  the  preceding  expression,  the  heat  penetration  depth  A1  is 
considered  a  function  of  time  only  and  so  is  A,,.  From  Equations 
1 1-3  and  1 1-7 ,  the  temperature  difference  is  given  as: 

Ti_co  -  Ti 3  (Vo.  ■  t2-»h1  '  F<VV1/n  +  (II"U) 
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Equation  11-12  then  becomes: 


(30,78)  g  _ 

(l  +  (k1/k2)(A2/A1  H  '§etTl-«  '  T2-»] 

a  T  3V8  1 

+  (T1-«  -  T2-co%  [i  +  (k1/k2)(A2/A1)J 
3A,/8  g2 

=  al  [TT  ■(■ kTkpTATS,  )']^T  tTl-»  -  T2-J 
3/(2A1) 

+  al[l  +  ( k, /k^)  U2/A, )  ]  [T1  -«  •  T2-»] 


With  Ti_oo  and  T2-«  9iven  by  Equations  II-l  and  II-2,  the  heat- 
balance  Equation  11-15  is  integrated  once  again  from  x  ■  o  to  x  * 
produci ng : 


3A,/8 

[1  +  (k1/k2)(A2/A1)] 


76 


3  r  3a^/8 

+  ^  |_i  +— (Tci/*^T(a27a^  ) 


3A,/  8 

“  “1[1  +  (k,/k2)(A2/A,)] 


(Tci^  -  75^) 
7a, 


3/(2A, ) 

*  “m  +  7^^4711 


(11-16) 
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There  is  another  equation,  similar  to  Equation  11-16,  obtained  by 
interchanging  the  indexes  1  and  2  in  Equation  11-16. 


To  solve  for  Aj  and  ^  as  t0  ^ow  they  change  with  time  6,  let: 

A]  =  S1^T[0  (11-17) 

A2  =  Sg/oigS  (11-18) 


where  and  S2  are  coefficients  to  be  determined, 
the  preceding  relations  is  the  fact  that  (A2/A-|)  is 
independent  of  time.  The  numerical  coefficients 
be  immediately  obtained  as: 


Implicit  in 
a  constant, 
and  Sg  can  then 


8/ [2  +  /agVouj  ] 


(II-19) 


S*  *  8/  [2  + 


(11-20) 


which  define  the  extents  of  the  two  boundary  zones  on  either  side  of 
the  interface.  At  the  interface,  the  temperature  is  given  by 
Equation  I I- 7  which  can  be  written  as: 


PVJ/t1 


+  P] 


(II-21) 


where: 


P  -  (kg/k^^/Ag) 

_ .  _  _ ,1/2 

=  (kg/k^ )*4.|/a2  {(2  +  /ag/ot^  )/(2  +  /ct^/otg)  J  (11-22) 


The  ratio  of  these  two  penetration  depths  (A^/Ag)  is  given  by: 


/a-j/a2 


2  +  /a^/a i 


r.1/2 


L  2  +  /a1 /a2. 


(11-23) 


1 1.3  THE  COMPLETE  TEMPERATURE  REPRESENTATION 


The  temperature  distributions  in  the  two  quarter-infinite 
regions,  1  and  2,  are  now  defined  by: 

(i)  Equations  II-l  and  I 1-2  remote  from  the  interface.  They 
are  of  course  the  well-known  error  functions  in  terms 
of  the  similarity  variables  x/2/a^S  and  x/2/5^5’. 

(1i)  Equations  II-3  and  II-4  near  the  interface  within  finite 
distances  A^  and  A 2  from  the  interface.  The  functional 
form  F(y1/A1)  or  F(y2/A2)  is  that  of  Equation  11-13,  and 
the  time-dependent  zone  depths  A-|  and  A2  are  given  by 
Equations  11-17  through  11-20. 

(Ill)  Equation  11-21  for  the  interface  temperature  T^.  The 
parameter  P  is  indicated  in  Equation  11-22. 
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In  order  to  bring  all  equations  to  a  common  basis,  a  similarity 
variable  n,  is  defined  as: 


^  2/a^9 


and  the  fully-established  temperatures  Tlco  and  are  expressed  by: 


Tl-«  a  er*cir\}) 


(11-24) 


T 


2-<» 


erfc 


(11-25) 


With  the  above  definitions  and  unifications,  numerical  results  are 
illustrated  for  the  following  combinations  of  the  physical  properties: 

*Vk2  “  5’  02^al  =  0*3 

for  which  the  parameter  P  becomes  0.298,  and  (A-|/A2)  *  1.48. 

In  Figure  II. 3  are  shown  the  established  temperature  distributions 
away  from  the  interface  and  the  interface  temperature  calculated 
from  Equation  11-21  for  the  parametric  values  specified  above.  The 
interface  zones  have  their  respective  temperature  variations  indicated 
in  Figure  II. 4  with  the  ordinate  values  i.e.  T^,  T^,  and  T2_w  chosen 
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'Z.  > 


Constant  Surface  Temperature); 


1.0  0.8  0.6  0.4  0.2  0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6 


Figure  I I. 4.  Cross-sectional  Temperature  Distribution  in  Interface 
Region.  Ordinate  values  for  ttj  *  .5  ( k-|/k2  *  5, 
a£/a1  «  0.3,  Constant  Surface  Temperature) 
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for  X/2.^"e  *  0.5,  while  the  abscissa  is  universally  valid  for  the 
parameters  so  chosen.  The  discontinuity  of  the  temperature  gradient 
at  the  interface  y  =  0  is  clearly  defined. 


I I. 4  CONSTANT  HEAT-FLUX  SOLUTION 


For  the  case  in  which  a  constant  (equal  in  both  regions)  heat 
flux  is  imparted  to  the  exposed  surface  x  =  0,  the  solution  method 
is  quite  similar,  but  with  the  exception  that  the  fully-established 
solutions,  instead  of  Equations  I I- 1  and  11-2,  are  now  given  by: 


(11-26) 


(11-27) 


where  Q  denotes  the  surface  heat  flux,  which  in  this  illustrative 
case  is  the  same  for  both  regions  (this  is  not  necessary  for  the 
method  to  work). 

Subsequent  developments  are  identical  up  to  and  including 
Equation  11-15  which  upon  integration  from  x  =  0  to  x  ■  «  yields 
the  following  equation: 
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(3A1/16)[(ot1/k1)  -  (a2/k2)]/n  +  (k^XA^)] 

+  *^§*(3/1 6) I (ot-|/k.| )  -  (og/kgJl^-tA^n  +  (k1/k2)(A2/A1 )]} 

*  a1(3A1/16)l(l/k1)  -  (l/k2)]/[l  +  (^/kgJCAg/A,)] 

+  a1(3/4A1)0[(a1/k1)  -  (a2/k2)]/[l  +  (k1/k2)(A2/A1 )]  (11-28) 

There  is  another  equation  obtained  by  starting  with  Equation  11-11  in 
region  2.  Integrating  Equation  11-11  from  y2  =  0  to  y2  =  «>,  an 
equation  similar  to  Equation  11-15  is  obtained.  Further  integration 
from  x  *  0  to  x  *  »  results  in  an  equation  which  is  a  companion  to 
Equation  11-28  and  can  be  obtained  from  Equation  11-28  by  interchanging 
the  indexes  1  and  2. 

It  is  of  importance  to  note  that  inter-region  conduction  is 
manifested  in  the  last  term  of  equation  11-15  in  the  form  of 
(T1<x>  -  T2  oo).  The  inter-region  conduction  influence,  which  modifies 
the  temperature  distribution  near  the  interface,  has  now  been  summed 
up  from  x  *  0  to  x  *  »  and  is  identified  in  Equation  11-28  as  being 
proportional  to  (a^/k^  -  a2/k2).  There  could  be  cases  depending 
on  the  thermo-physical  properties  of  the  media  so  that  the  integrated 
sum  of  the  Inter-facial  heat  conduction  becomes  zero.  This  occurs 
when  p^c1  *  pzc2’  F*9ure  11,5  illustrates  this  particular  combination 
which  is  due  to  a  reversal  of  the  transverse  heat  conduction.  At  or 
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near  this  condition,  the  analysis  of  axial  heat-balance  integral 
would  be  quite  inaccurate  and  a  modification  of  the  method  is 
necessary. 

The  cause  of  trouble  as  revealed  in  Figure  I I. 5  leads  to  an 
obvious  remedy.  Instead  of  integrating  Equation  11-15  and  its 
companion  (obtained  by  interchanging  the  indexes  1  and  2  in  Equation 
11-15)  from  x  *  0  to  x  =  »,  two  heat-balance  integrals  are  obtained 
from  each  of  these  two  equations.  One  is  obtained  from  x  *  0  to 
x  =  Xq,  the  cross-over  point  in  Figure  I I. 5  and  the  other  is  obtained 
from  x  *  xQ  to  x  *  «>.  In  this  way,  cancellation  of  the  transverse 
conductance  is  prevented,  and  the  conduction  effect  is  duly  accounted 
for  in  each  range.  In  each  range  then,  variations  of  and  a0 
with  time  are  determined.  It  is  to  be  noted  that  the  time-variation 
of  A-j  or  A2  in  the  first  range  (x  =  0  to  x  ■  xQ)  may  be  at  variance 
with  the  time-variation  in  the  second  range  (x  >  x^).  The  discrepancy 
between  the  values  on  either  side  of  the  cross-over  point  is  not 
detrimental  to  a  smooth  transition  of  the  temperature  at  the  cross¬ 
over  point,  since  at  the  latter  there  is  no  inter-region  temperature 
difference  or  in  other  words  the  heat  penetration  depths  A-j  and  A2 
are  moot  at  this  location. 
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III.  PARALLEL  HEAT  FLOW  IN  LAMINATED  COMPOSITES 


Consider  a  laminated  composite  consisting  of  two  alternate  layers 
of  two  different  materials  and  occupying  a  half-space  x  >  0.  For  heat 
flow  along  the  x-direction,  a  unit-cell  can  be  isolated  from  the 
assembly  and  the  cell  consists  of  a  half-lamination  of  material  1 
adjoint  to  a  half-lamination  of  material  2.  Figure  III.l  depicts 
such  an  arrangement. 

The  boundary  conditions  on  the  symmetry  surfaces,  y  *  0  and  y  =  b, 
are  of  course  zero  heat-conduction.  The  interface  boundary  conditions 
are  equal  temperature  and  equal  heat  flux.  The  boundary  condition 
at  the  exposed  surface  x  =  0  will  be  taken  as  (i)  that  of  a  constant 
heat  flux,  or  ( i i )  that  of  a  constant  surface  temperature.  The  composite 
body  is  initially  at  zero  temperature. 

Before  proceeding  with  various  approximate  methods,  exact  analytical 
solutions  were  first  obtained  for  laminated  composites  with  a  finite 
length  in  the  principal  heat  flow  direction.  Instead  of  extending  to 
infinity  along  the  x-direction  in  Figure  II. 1,  the  composite  terminates 
at  x  *  L. 

Two  boundary  conditions  at  x  a  L  were  used  alternately  in  that 
analytical  study,  namely  (i)  a  perfect  insulation  and  (ii)  constant 
temperature  (equal  to  its  initial  value).  Hence,  all  together  four 
analytical  solutions  expressed  by  infinite  series  were  developed  and 
presented  in  Appendix  C. 
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Insulated 


Figure  III.l.  Schematic  of  a  Laminated 
Composite 
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The  reasons  for  that  effort  are  two-fold.  First,  exact  solutions 
for  composite  configurations  of  practical  interest  are  sparse,  and 
secondly  the  availability  of  solutions  for  the  geometry  described 
would  establish  a  standard  with  which  approximate  methods  can  be 
compared  from  the  viewpoint  of  labor- vs. -accuracy  trade-off.  Even 
though  no  numerical  results  have  come  forth  from  that  effort,  because 
of  the  resulting  complexity,  they  are  presented  in  Appendix  C  for 
completeness  and  possible  future  uses. 

Turning  now  to  the  approximate  methods  of  analysis,  the  geometry 
depicted  in  Figure  1 1 1 . 1  is  taken  up  next.  Conceptually,  the  temperature 
variations  in  both  regions  are  expected  to  be,  at  small  times,  very  much 
like  those  in  a  semi-infinite  domain  occupied  by  each  constituent  with 
the  other  absent  in  the  half-space.  As  an  example,  the  temperature 
responses  for  the  case  of  a  step-rise  in  their  surface  temperature  at 
x  *  0  are  given  by: 

T-j  3  erfc(x/2^o|e') 

Tg  *  erfc(x/2v6^T) 
for  small  times. 

As  time  proceeds,  inter-region  heat  transfer  takes  place  which  is 
manifested  by  transverse  temperature  variations  which  result  in  a 
significant  deviation  from  the  solutions  just  described. 

The  preceding  description  of  how  the  temperature  va-ies  qualitatively 
can  be  replicated  by  adopting  the  concept  that  in  each  region  the 
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temperature  distribution  is  comprised  of  two  parts:  a  first  part  which 
is  dependent  on  the  axial  distance  x  and  time  6  only,  and  a  second  part 
which  is  a  product  of  two  functions,  one  dependent  on  (x,  8)  and  the 
other  on  (y,  0). 

Obviously  at  the  interface  y  *  a  (Figure  II . 1 ) ,  these  two  regions 
must  have  a  common  temperature  T...  Such  a  temperature  picture  is  given 
by  the  following: 

T,  -  t,  +  (t.  -  t1)F(y,  9)  (y  <  a) 

T2  =  t2  +  (t..  -  t2H(y,  6)  (a  <  y  <  b) 

where  t| ,  1 2,  ti  are  functions  of  (x,  6). 

Two  crucial  roles  are  to  be  fulfilled  by  F(y,  9)  and  $(y,  9). 
First,  at  the  interface  y  =  a,  these  two  functions  must  assume  a 
value  of  unity  so  that  Tj  *  T2  *  t.. .  Secondly,  at  time  zero,  both 
functions  must  vanish.  These  two  requirements  produce  a  singularity 
at  0  a  0  and  y  *  a,  and  are  therefore  met  by  the  following  functional 
form  for  F: 

F(y,  0)  *  cosh8-|y/coshe-|a 

where  8j  Is  an  Implicit  function  of  time  0,  such  that  the  singular 
behavior  of  8^  ■»  ®  at  9  ■*  0  yields  the  desired  characteristics  of 
F(y,  0)  as  discussed  before.  Of  course,  the  fact  that  cosh  8^y  has  a 
zero  slope  at  y  ■  0,  the  symmetry  line,  is  a  built-in  feature. 
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Similarly,  the  function  4>(y»  9)  Is  given  the  form: 

<j>(y,  9)  3  coshS2(y  -  b)/cosh62(a  *  b) 

which  has  a  zero  slope  at  y  =  b,  thus  fulfilling  the  insulated  condition 
at  y  *  b. 


III.l  FUNCTIONAL  REPRESENTATION  OF  THE  TRANSVERSE  TEMPERATURE  VARIATION 


Consider  the  following  functions  to  denote  the  temperatures  in 
regions  1  and  2  respectively. 


,  cosh8,y 

T1  =  h  +  "  tl)cosh81a 

for  region  ?,  and 

,  cosh8?(y  -  b) 

T2  3  t2  +  TT T(tl  "  Vcoshe2(a” b) 

for  region  2. 


(III-D 


( I I I- 2 ) 


The  parameter  A  is  defined  by: 

/k2\  / S2\  tanh02(b  -  a) 

X  *  \^l" )  \®i" )  tanK^a 


(II 1-3 ) 


The  second  terms  of  Equations  III-l  and  I I 1-2  represent  the  interaction 
of  one  region  with  another  through  conduction  across  tne  interface  y  s  a. 
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If  the  second  terms  are  zero,  then  tj  or  tg  will  be  the  temperature 
distribution  in  each  region  when  the  other  is  absent.  These  temperature 
distributions  in  such  an  extreme  case  are  of  course  known  from  previous 
work  and  can  therefore  be  used  in  the  construction  of  these  one¬ 
dimensional  solutions  t^  and  tg. 

In  these  equations,  8^  and  3g  are  functions  of  time  and. can  be 
thought  of  roughly  as  1/^0  and  l/v5^e  respectively,  although  not 
precisely  equal  to  them.  Thus,  when  at  small-time  instants,  81  and  8g 
tend  to  infinity,  the  second  terms  of  Equations  III-l  and  I I 1-2  have 
a  vanishingly  small  influence.  In  this  way,  the  temperature  responses 
in  regions  1  and  2  independent  of  the  other  are  preserved. 

The  interface  conditions  at  y  =  a  are  fulfilled  by  their 
functional  formations.  These  conditions  are: 


kj(3Tj/3y)  =  kg(3Tg/3y) 

and  are  satisfied  by  Equations  II-l,  I 1-2  and  III-3.  Furthermore,  the 
insulated  condition  at  y  3  a  and  y  a  b  are  also  complied  with. 

In  summary.  It  should  be  noted  that  t1  and  tg  are  functions  of 
x  and  6,  and  that  T^  and  Tg  give  the  spatial-time  distributions  of 
temperature  in  these  two  regions.  Undetermined,  however,  are  the 
time-parameters  8^  and  8g  as  well  as  t^  and  tg.  It  is  the  purpose  of 
the  heat-balance  Integral -approach  to  obtain  these  variations. 
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1 1 1. 2  THE  HEAT-BALANCE  EQUATIONS 


8y  averaging  the  temperature  in  each  region  as  given  by 
Equations  III-l  and  III-2,  there  result: 


(II 1-4) 


( 1 1 1-5 ) 


The  equations  which  govern  and  Tg  are  obtained  by  averaging  the 
governing  differential  Equations  I 1-6  and  II-7.  In  carrying  out  the 
process  of  averaging,  the  Insulated  (zero  conduction)  condition  at  y  *  0 
and  y  ■  b  and  the  equal  heat  flux  condition  at  y  *  a  are  utilized.  The 
governing  equations  after  being  averaged  then  become: 


3t,  a  T  i  tanhS^a  T  , 

8l*  J  +  tt 


tanh&,a 


tanh6,a  . 

4(t,  -  t,] 


X  B^a  30l,*2  1 


a2  r  .  tanh8,a- 

a!  +  (t2  ’  V  TTJ  6,3 


f i  ~r  ,  (t-  -  t,  )B.,tanh(3ia"| 


(Iii-e: 


h2  .  *  x  3  r  i  tanh®2(b " 

36  ‘  V  38  [TTT  g2tBT“a) 


^  tanhe2(b  -  a)  3 
+  TTT  e2(b  -  a)  -55-  [t1 
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(ti  * 


U 


tanh62(b 


a) 


2'  1  +  X  62(b  -  a) 


TF^aT 


■ 

1 

nr 

i» 


(ti  “ 


t2)e2tanh02(b  -  a) 


(III-7) 


In  the  preceding  equations,  the  second  terms  on  the  left-hand 
sides  denote  the  sensible  heat  as  the  temperature  profiles  in  the 
transverse  direction  (i.e.,  the  y-direction)  change  with  time.  The 
reason  for  this  identification  is  that  the  second  terms  calculate 
the  time-change  of  ^  and  @2  respectively.  Thus,  the  three  terms 
on  the  left-hand  sides  can  be  interpreted  as  follows:  the  first 
terms  represent  sensible  heat  change  due  to  a  bulk  temperature  change; 
and  the  third  terms  represent  the  change,  not  of  the  profiles  as  do 
the  second  terms,  but  of  the  bulk  value  (t-,  -  t0). 

On  the  right-hand  sides  of  these  equations,  it  is  easy  to 
recognize  that  the  first  terms  reoresent  the  net  axial  conduction  and 
the  second  terms,  the  transverse  or  Interface  conduction. 

In  laminated  composites,  the  transverse  dimension  a  or  b  Is 
usually  much  smaller  than  the  axial  dimension  -  the  direction  which 

coincides  with  the  principal  heat  flow  path.  Thus,  a  diffusion  time 

2  2 

for  the  transverse  direction  is  naturally  b  or  b  /a2,  which 
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characterizes  the  time  scale  the  transverse  temperature-change 
takes  place.  There  is,  of  course,  a  diffusion  time  for  the 


2  2 

longitudinal  direction  indicated  by  z /a^  or  z  /ot2,  which  is  indicative 


of  the  time  scale  for  heat  propagation  in  that  direction.  In  the 
region  where  x  »  b,  the  temperature  change  (with  time)  is  much 
slower  in  the  x-direction  than  that  in  the  transverse  direction,  due 
to  Interface  conduction  across  a  much  shorter  distance.  From  this 
dimensional  argument,  it  is  plausible  to  equate  the  second  term  on 
the  left  hand  side  of  Equation  I I 1-6  with  the  second  term  on  the 
right-hand  side  of  it.  Of  course,  the  same  can  be  said  of  Equation 
1 1 1-7.  By  such  isolation  from  the  overall  equations,  the  "transverse- 
balance"  equations  are  now: 


»ici  I?  [  1  h  ta?a]°  r-  [rhr  <ht*nhv] 


(II 1-8 


P2C2 


j  tanhB2(b  -  a)" 

k2 

1  +  X  82(b  -  a) 

(b  -  a) 

^  ^  P2wmnP2\u  •  »/ 

( 1 1 1 -9 


Because  of  Equation  1 1 1-3,  which  equates  the  Interface  heat  conduction, 
the  right-hand  sides  of  Equations  I I 1-8  and  III— 9  are  equal  (except 
for  the  factors  a  and  b  -  a).  A  few  manipulations  show  that 
although  8^  and  82  are  both  unknowns  at  this  stage,  they  are  related 
by  the  following  simple  relation: 


*  •/0|/a g 


(III-10) 


and  Equation  I I 1-3  for  X  can  be  re-cast  as: 

_ tanK?,{[(b/a)  -  1  ]/a,/a5} 

x  ‘  ^  »»-"> 
In  both  expressions: 

8^  3  6^,  and  f2  =  (III-12) 

The  above  developments  serve  the  purpose  that  only  one  of  the  two 
equations  Involving  8-,  and  $2  needs  to  be  solved.  Thus,  Equation  III-8 
will  be  sufficient  to  determine  8-|  or  ^  as  a  function  of  time  e. 

The  pertinent  boundary  condition  Is  at  0  -*•  0,  8^  -*■  »  which  is  required 
in  the  temperature-profile  Equations  III-l  and  III-2.  For  8^  ®, 

conslstant  with  $2  +  the  second  terms  of  Equations  III-l  and  II 1-2 
both  vanish  as  they  should,  at  0  •+  0. 


III. 3  THE  TRANSVERSE  HEAT-BALANCE- INTEGRAL  EQUATIONS 

Regrouping  the  pertinent  equations  Into  this  section,  and  using 
the  following  simplified  notations: 
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01  * 


e2  =  s2a 


0  *  (o^e/a^). 


(III-12) 


the  equation  to  be  solved  for  ^  Is, 


d  r^_  tanl,ei  1 .  _x_ ,  tanh5 

dff  V  *  »  ?,  J  '  1  +  *  6'  6' 


(II 1—8 ) 


with  the  boundary  condition  0  -*•  0,  g^  -►  ».  The  parameter  A  is 
defined  by: 


•4 


X  »J(k2P2c2/k1P1c1)  tanh  J  0^—  (- -  1) 


t  n  /b 


)  j  /tanhgj 


(III-ll) 


Defining  a  new  symbol  H  such  that: 


tanhg, 

LI  —  A  * 

H  -TTT-J— 


(III-13) 


equation  I I 1-8  can  be  put  as: 


4*'  He* 
de  1 


(III— 14) 


35 


with  the  boundary  condition  of  8-|  ■  °»  at  0  ■  0.  Appendix  A  shows 
the  numerical  method  developed  to  solve  Equation  III-14  as  well  as 
the  asymptotic  continuation  at  large  values  of  'e. 

For  kg/kj  ■  10,  ctg/oj  *  3  and  b/a  =  2,  the  variation  of  ^  with 
"9  Is  shown  in  Figure  III. 2.  At  small  values  of  9",  frj  =  1  Hence 

for  small  times,  the  transverse  temperature  variation  as  shown  in 
Equation  1 1 1-1  has  the  form  of: 

cos  hey  a  CQSh^i  (y/a)  _  COsh(v/a)/^f 
C°Sh6a  cosh^l  cosh  (1/^t) 

for  the  laminate  in  region  1,  i.e.,  y  <  a.  At  a  short  distance  from 
the  interface  y  =  a,  let  y‘  =  a  -  y,  the  transverse  temperature 
variation  as  given  above  can  be  expressed  by  the  following  (for 
-  «): 

cos h3y  .  -6(a  -  y)  _  -6y'  .  -(y'/a)/^! 
coshea  e  e  e 

which  shows  the  correct  combination  of  the  similarity  variable  in  heat 
conduction  theory. 
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1 1 1. 4  THE  AXIAL  HEAT  DIFFUSION  EQUATIONS 


In  order  to  solve  Equations  I I 1-6  and  III-7,  which  are  the 
results  of  averaging  the  diffusion  equations  along  the  transverse 
direction  over  these  two  regions,  the  averaged  Equations  I I 1-6  and 
1 1 1-7  can  be  expressed,  with  7  *  x/a,  as: 


( III-15) 


(III-16) 


where  H  Is  defined  by  Equation  I I 1-13.  The  definition  for  R  in  the 
preceding  equation  is: 


(pc), 

3  (pc)2(J  -  1) 


(III-17) 


which  amounts  to  the  ratio  of  the  thermal  capacities  In  regions 
1  and  2. 

In  these  two  equations,  the  coefficients  of  the  second  terms  on 
the  right-hand  side  are  of  course  presumed  known,  as  they  have  been 
obtained  from  Equation  II 1-14.  Here,  the  focus  Is  on  the  average 
temperatures  T,  and  Tg  which  are  the  results  of  axial  diffusion  and 
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transverse  Inter-region  heat  conduction.  The  merit  of  only  having  to 
solve  Equations  III-15  and  III-16,  in  contrast  to  solving  the  full 
diffusion  equations  for  these  two  regions,  lies  in  the  fact  that  the 
transverse  coordinate  y  disappears  from  these  two  reduced  equations 
by  virtue  of  the  transverse  time-dependent  parameters  6^  and 

This  seemingly  minor  simplification  of  the  equations  results  in 
a  very  significant  reduction  In  labor  required  to  obtain  the 
temperature  distributions.  The  key  lies  in  a  proper  choice  of  the 
profile  functions  for  F(y,  0)  and  <p(y,  0)  in  these  two  regions. 

Perturbation  approaches.  In  order  to  address  the  various 
solution  techniques  available,  let  the  inter-region  conduction 
coefficient  be  expressed  by  the  following: 

H6? 

L  ■  t  -.-imi  ♦ ») 

and  Equations  I I 1-1 5  and  III-16  are  re-shaped  to: 

aT.  a2T, 

— !■  *  — 4-.  l(T,  -  T.)  (HI-19) 

aF  ax'1 

—  •(?)—§  -  RL(T2  -  T})  ( 1 11-20) 

aF  W  aT2  2  1 
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The  Inter-region  conductance  L  has  the  following  characteristics,  as 
deduced  from  asymptotic  analyses  (see  Appendix  A): 

L  *  lA/f  for  TF  -»■  0 

L  «  constant  f or  "5  -  « 

Also  of  significance  is  the  thermal  capacity  ratio  R,  defined  by 
Equation  I I 1-17. 

Conceptually  then,  the  system  of  Equations  III-19  and  I 11-20 
with  proper  boundary  conditions  at  x  *  0  can  be  solved  by  various 
perturbation  techniques  (the  use  of  the  Laplace  transforms  does  not 
obviate  the  use  of  perturbation,  as  the  transform  parameter  is  not 
simply  obtained).  Therefore,  one  approach  is  to  express  the  solutions 
and  T2  in  series  with  time  IT  as  the  perturbation  parameter,  and  the 
scheme  is  to  determine  the  terms  in  these  two  series  on  a  term-wise 
basis.  Thus  by  writing  T.  and  Tg  as: 

Tj  «  Gq(x,  0)  +  V^Cx.  <T)  +  0G2(x,  0")  +  .  .  .  . 

T2  ■  gQ(x,  0)  +  Vfg^x,  IT)  +  0g2(x,  ?)  +  .... 

where  the  first  term  solutions,  G0  and  g0,  are  the  infinite-region 

solutions  without  inter-region  conduction.  This  is  true,  of  course, 
for  problems  with  initial  temperature  zero  in  both  regions. 
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A  second  perturbation  approach  Is  through  the  thermal  capacity 
ratio  R.  If  R  is  very  small,  which  In  a  limiting  sense  may  be  put 
to  zero,  then  Equation  I 11-20  can  be  solved  immediately.  In  other 
words,  the  temperature  distribution  T 2  will  be  unaffected  by  Tj.  Using 
this  solution  of  T2,  Equation  III-19  can  be  solved  for  T^,  which  is 
used,  in  turn,  to  solve  for  T2  from  Equation  III-18.  Such  a  procedure 
is  equivalent  to  an  iteration  approach.  Alternately  of  course,  a 
formal  series  expansion  in  terms  of  R  can  be  used,  such  as: 

Tj  s  SqCx.IT)  +  RS-jCx,  IF)  +  .  .  . 

T2  =  i|»  (x,  e-)  +  R«p1  (x ,  0)  +  .  .  . 


Double  Heat-Balance  Integral -Method.  Equations  III-19  and  I 11-20 


are  the  results  of  using  heat-balance  integrals  in  the  transverse 
direction;  they  are  considerably  simpler  than  the  full  diffusion 
equations.  Under  certain  conditions,  these  two  equations  can  be 
approximately  solved  by  using  once  more  heat-balance  integrals  in 
the  axial  direction.  The  resulting  approximate  solutions  for  the 
two  regions  may  be  termed  the  double  heat-balance  integral -method. 

To  demonstrate  its  efficacy  and  deficiency,  two  types  of  surface 
heating  will  be  considered:  constant  surface  temperature  rise  at 
x  ■  0,  and  constant  surface  heat  flux.  For  the  former  type  of 
boundary  condition,  the  method  of  double  heat-balance  integrals  shows 
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a  significant  advantage,  while  for  the  latter  type  of  boundary 
condition,  some  limitation  of  the  results  must  be  observed. 


Constant-Temperature  Case.  First,  consider  the  problem  of  a 
constant  surface  temperature  rise.  Let  the  surface  temperature  at 
x  »  0  of  both  regions  be  impulsively  raised  to  a  finite  value  from 
their  initial  temperature  of  zero.  If  there  is  no  inter-region 
conductance,  i.e.,  L  =  0,  then  the  solution  for  each  region  is 
independent  of  the  other  and  is  given  by  the  following  for  region  l: 

=  erfc(x/2Vt)  (III-21) 

For  region  2,  the  solution  is: 

(111-22) 

Essential  to  the  method  of  heat-balance  integral  is  a  proper 
choice  of  temperature  profiles  for  these  two  regions  when  inter- 
region  conduction  is  to  be  accounted  for.  Accordingly  then,  the 
error-functions  are  chosen  to  represent  the  temperature  profiles. 

But  in  liru  of  two  time-dependent  parameters  Yj  and  Yg  are 
embedded  in  the  preceding  expressions;  and  the  chosen  profiles 
are  now: 
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a  erfc(y^x) 


(III-23) 


T2  *  erfc(y2x) 

in  which  the  variations  of  y^  and  y2  versus  time  O’  are  to  be 
determined  via  the  heat-balance  integral  method.  At  very  small  times 
IT,  Equations  1 11-23  and  1 1 1-24  must  replicate  the  independent  solutions. 
Equations  111-21  and  III-22,  as  the  inter-region  conductance  is 
initially  absent.  This  requirement  establishes  the  necessary 
boundary  conditions  on  y-j  and  y2,  i.e.  for  J  -*■  0: 

y1  - 

y2»^7^/(2Vt)  Cl  11—24 ) 


In  fact,  a  more  appropriate  description  is: 

Yj  »  •  at  o'  -*•  0 

Y2  3  «  at  ¥  -*■  0 

and  their  singular  behaviors  at  8  -*■  0  are  determined  by  the  heat- 
balance  integrals. 


(II 1-25 ) 
(II 1-26 ) 
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The  next  step  is  to  integrate  Equations  I I 1-19  and  I I 1-20  from 
x  *  0  to  x  *  »  with  the  temperature  profiles  Tj  and  T2  given  by  Equations 
111-23  and  II 1-24.  The  result  is  two  equations  with  and  Y2  as  the 
new  dependent  variables. 


(II 1-27) 


(II 1-28 ) 


The  associated  boundary  conditions  are  those  prescribed  by  Equations 
1 1 1-25  and  1 1 1-26.  It  is  to  be  recalled  that  the  time-dependent  transverse 
conductance,  L,  having  been  determined  by  virtue  of  solving  Equation  III-13, 
is  a  Dre-determined  function  in  Equations  I I 1-27  and  I 11-28. 

For  starting  the  solution  of  these  two  equations,  the  initial 
boundary  conditions  of  Equations  I 11-25  and  I 11-26  can  be  satisfied 
by  taking: 

—  *  +  •  .  . 

Y1 
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The  numerical  coefficients  s1  and  s2  are  immediately  shown  to  be 
s1  a  2  and  s2  =  which  render  the  asymptotic  solutions  (for 

IF  -*•  0)  the  same  as  the  independent  solutions.  As  time  goes  on,  the 
influence  of  the  transverse  conductance  L  starts  to  modify  the 
solutions.  Implementation  of  the  solution  can  proceed  by  a  simple 
integration  method.  Results  obtained  in  this  manner  are  discussed 
in  Section  V  in  comparison  with  other  solution  techniques. 


Constant  Heat-Flux  Case.  In  this  case,  the  temperature  profiles 
in  the  two  regions  are  guided  by  their  respective  independent  solutions. 
Thus,  starting  with  the  profiles: 

T,  -  [aQ/k1y1  ]  ierfc (7^  (III-29) 


T2  =*  [aQ/k2Y2]  ierfc(y2x) 


(II 1-30) 


where  y-j  and  y2  are,  as  in  the  case  of  constant  surface 
the  time-dependent  parameters. 

Integrating  Equations  I I 1-19  and  I I 1-20  from  7  »  0 
with  7^  and  T2  as  the  profiles,  the  resulting  equations 


temperature, 

to  x  =  ~ 
become: 


(II 1-31 ) 
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(II 1-32) 


The  preceding  equations  can  be  solved  in  a  closed  form  for  the 
boundary  conditions  Yi  ♦  ••  ?2  *  00  at  ^  "** 


(III- 33 ) 


which,  when  introduced  into  Equation  III-31  or  III-32,  gives  direct 
solutions  for  y^  or  The  surface  temperature  variations  (with 
time)  for  these  two  materials  are  given  by: 

T,  -  IaQ/k^j/7,  (II 1-34) 

T2  «  taQ/kj/jrHkj/kgl/yg  (HI-35) 

Numerical  results  obtained  from  the  preceding  equations  are  discussed 
in  Section  V  In  conjunction  with  results  obtained  by  other  methods. 
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To  understand  the  nature  of  the  double  heat-balance  integral 
method.  It  is  necessary  to  examine  Equations  III-19  and  1 11-20,  which 
contain  the  inter-region  conductance  effect.  The  inter-region 
conductance  exhibited  in  these  equations  is  of  a  local  character, 
i.e.,  at  an  axial  location,  x.  When  axial  heat-balance  integrals  are 
applied,  the  inter-region  conductance  is  therefore  summed  up.  Thus, 
the  local  variation  of  L  is  thereby  lost.  This  is  a  source  of 
inaccuracy  in  the  final  results.  In  an  extreme  case,  the  inter- 
region  heat  conduction  may  reverse  itself,  i.e.,  from  x  s  0  to  x  *  x, 
(T2  -  )  >0  and  from  x  *  x^  to  x  «  »,  (T2  -  T^)  <  0.  The  average 

of  the  interfacial  conduction  may  be  quite  influential  on  the  final 
results. 


Finite-Difference  Solution.  The  method  of  perturbation  described 
in  Section  III. 5.1  is  theoretically  exact  and  conceptually  elegant. 
From  a  practical  point  of  view,  the  method  Is  tedious,  for  it 
requires  a  number  of  terms  in  order  to  achieve  some  degree  of 
accuracy.  Similarly,  the  use  of  the  Laplace  transform  -  again 
conceptually  exact  -  does  not  render  Itself  to  immediate  numerical 
extraction,  because  the  transform  parameter  does  not  have  simple 
roots.  Consequently,  the  use  of  small-time  expansion  or  large-time 
expansion  in  the  inversion  procedure,  which  constitutes  a  de  facto 
perturbation  process,  negates  the  conceptual  advantages  the  Lapl ace- 
transform  method  possesses. 
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On  the  other  hand,  the  double  heat-balance  method  presented  In 
Section  III. 5. 2  reduces  the  analytical  labor  to  a  minimum,  but  In  so 
doing  the  local  (along  the  axial  position)  Inter-region  conductance 
effect  Is  averaged  out  or  may  be  mutually  cancelled.  Only  the  total 
Inter-region  transverse  conduction  Is  recognized,  as  Is  pointed  out 
previously  In  Section  II. 

Between  these  two  approaches  lies  the  method  of  solving  the 
axial  diffusion  Equations  III-19  and  I 11-20  by  a  finite-difference 
procedure,  which  in  a  broader  sense  belongs  to  the  exact-solution 
category.  The  advantage  of  solving  the  axial  diffusion  equations 
over  the  full  diffusion  equations  in  two  spatial  coordinates  is 
readily  apparent  if  it  is  recognized  that  the  labor  of  the  former 
is  about  one-tenth  of  the  latter.  For  this  reason,  the  axial 
diffusion  Equations  III-19  and  III-20  as  well  as  the  corresponding 
equations  for  fibrous  composites  are  used  extensively  in  this  report. 

Appendix  B  contains  a  synopsis  of  the  finite-difference 
formulation,  and  the  results  therefrom  are  presented  in  Section  V 
together  with  those  from  other  methods. 
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IV.  PARALLEL  HEAT  FLOW  ALONG  UN I -DIRECTIONAL  FIBERS 


Consider  an  idealized  dispersion  pattern  of  uni-directional 
fibers*  in  a  matrix  material  such  that  the  fiber  centers  form  an 
equilateral  triangular  pattern  as  depicted  in  Figure  IV. 1,  which 
shows  a  cross-sectional  view. of  such  a  composite  material.  When 
planar  heat  flow  takes  place  along  the  fiber  axial  direction,  it 
is  possible  to  visualize  a  repeating  cell  in  which  the  thermal 
response  is  representative  of  the  entire  body. 

Such  a  unit  cell  is  a  hexagon  with  a  fiber  located  in  its 
geometrical  center  as  indicated  in  Figure  IV. 1.  Intuitively,  the 
hexagonal  enclosure  of  the  fiber  can  be  replaced  by  a  circular 
enclosure  containing  the  same  amount  of  the  matrix  material ,  and 
the  boundary  condition  on  the  circular  enclosure  is  equivalent  to 
that  on  the  hexagonal  enclosure,  i.e.,  an  insulated  condition. 

Mathematically  speaking,  a  great  simplification  has  been 
achieved;  both  boundaries  are  now  expressed  in  circular  coordinates 
while  preserving  almost  complete  thermal  equivalence.  Hence,  the 
problem  is  transformed  into  two  circular  cylinders  with  the  matrix 
material  enclosing  the  fiber  as  is  now  illustrated  in  Figure  IV. 2. 


*Here  and  elsewhere  In  this  report,  the  terms  fiber  and  tow  are 
used  interchangeably. 


Equivalent  Circle 


Figure  IV. 2.  Equivalent  Two-cylinder  Configuration 
for  a  Uni-directional  Fiber-composite 
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To  preserve  parallelism  with  the  case  of  laminated  composites* 
the  radii  are  named  a  and  b  and  the  principal  heat  flow  Is  along  the 
axial  direction  x.  For  Illustrating  the  heat-balance  integral 
approach,  consideration  Is  given  to  the  problem  of  a  uniform  Initial 
temperature  (zero)  and  an  impulsively  applied  heat  flux  Q  on  the 
surface,  x  »  o.  The  material  properties  are  ,  p.j,  c^for  the 
fiber  region  and  kg,  P2«  c^for  the  matrix  (annular)  region. 

With  the  simplified  configuration  shown  in  Figure  IV. 2,  the 
problem  Is  the  same  as  that  treated  In  Section  III  except  here  a 
circular  configuration  replaces  the  planar  configuration  in 
Section  III. 


IV. 1  FUNCTIONAL  REPRESENTATION  OF  THE  TRANSVERSE  VARIATIONS 


A  mathematical  function  which  reasonably  reflects  the  spacial 
distribution  of  the  temperature  must  contain  r-variation  and 
x-varlation.  Expressed  in  a  product-form,  the  temperatures  In  regions 
1  and  2  can  be  taken,  similar  to  Equations  III-l  and  III-2,  as: 


Ti  *  *i  +  (ti 


[coshS,r 

- '-1 

cosh^a  _ 


T-  ■  t,  +  (t.  -  tJ 


’cosh82(rz  -  bz) 
,2"*2  '  ^  ’  W2#  [cosh82(a2  -  b*). 


(IV-1) 


( IV— 2) 
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The  radial  variation  (In  space  and  in  time)  is  expressed  through 
the  terms  In  brackets.  It  is  readily  seen  that  the  conditions  of 
zero  temperature-gradient  at  r  *  0  and  r  =  b  are  fulfilled  by  the 


form  so  chosen.  Furthermore,  at  r  *  a,  the  interface,  Tj,  is  equal  to 
T2  as  is  required.  The  equal  heat-flux  conditions: 

k1(3T1/3r)a  =  k2(3T2/3r)a 

\ 

result  in  an  expression  for  T^ ,  the  interface  temperature,  as: 

tj  *  +  Xt2)/(1  +  X)  ( IV- 3 ) 


and 


k2\/e2\  tanh02(b2  -  a2) 
M'V  tanhB^a2 


(IV-4) 


which  is  parallel  to  that  for  the  planar  case,  see  Equation  III-3. 
The  profiles.  Equations  IV- T  and  IV-2,  can  be  re-cast  as: 


T 


1 


T 


2 


^1  *  ^2  " 


t2  +  (t'j  - 


V 


t2) 


2 

x  r  cosier  I 

1  +  x  Lcosh^a2  J 

1  rcosh62(r2  -  b2)"| 

1  +  X  Lcosh62(a2  -  b*)J 


(IV-la) 


(IV-2a) 
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Temperature  variations  in  x  and  8  are  embedded  in  t^  and  t2, 
whereas  temperature  variations  in  r  and  8  are  expressed  through  8-j , 

82,  r»  and  X.  The  time-dependent  parameters  6-j  and  82  determine 
the  transverse  profiles  in  both  materials. 

As  in  planar  case.  Section  III,  it  is  of  importance  to  note 
that  the  boundary  conditions  on  8-j  and  B2  are  that  as  8  -*•  0,  and 
B2  approach  infinity.  The  latter  condition  yields  a  singularity  at 
the  interface  r  »  a,  while  it  renders  the  temperature  in  the  transverse 
direction  uniform  in  their  respective  cross-sections. 


IV. 2  THE  HEAT-BALANCE  EQUATIONS 


First,  the  temperature  profile  Equations  IV-la  and  IV-2a  are 
averaged  transversely  in  each  region,  resulting  in: 


'2 


*1  +  (t2  -  t 


J  x  tanhe/j 


2  .2, 


r  «.  tanhS9(b  -  a  )  -i 


(IV-5) 


( IV— 6 ) 


Next,  the  governing  heat  diffusion  equations  are  averaged  in 
the  r-dlrection  for  each  region.  The  resulting  governing  equations 
are  in  the  following  form,  similar  to  those  in  Section  III: 


54 


it  1 1 3  r  *  tanhB?a2i  .  r * 
<2  L  rr- 


tanh6,a  -|a 

I^^-h(t2"tl) 


2T  r 

"=?  +  (t2  _  *1  }  [r+ T  (61a2)tanh(e]a2)  J 


( IV— 7 ) 


tanhe2(b2  -  a2) 
82(b2  -  a2) 


1  tanhe2(b2  -  a2 
~  ez(b2  -  a2) 


J]3 

(t*  •  V 


’-&*  lt'  '  t2)['  ^((b/a)2-l)tanhS2<b2  '  a2)J  (IV'8) 


( ( b/ a )  -  1) 


where  0  *  a^e/a  and  x  *  x/a. 


IV. 3  THE  TRANSVERSE  HEAT-BALANCE  INTEGRAL  EQUATIONS 


As  Is  explained  In  Section  III,  the  interface  conduction 
represented  by  the  second  term  on  the  right-hand  side  of  Equation 
IV-7  or  IV-8  is  associated  with  the  time-rate  change  of  the  transverse 


1 


temperature-profile  represented  by  the  second  term  on  the  left-hand 
side  of  the  equation.  These  two  terms  can  be  isolated  from  each 
equation.  Such  a  "transverse  equilibrium"  gives  the  following  heat- 
balance  integrals: 


4. r  >  uny  1 >  a 


2  2 
a  tanhB^a 


g,  d  r  ,  -1)1  ,  2  ?,  b  '‘ 

oj  jg-  L  1  +  *  J  '  4  TTT  B2a  tanhB2a  {  y 


(IV-9) 


(IV-10) 

Because  of  the  definition  of  X  in  Equation  IV-4,  the  right-hand  side 
terms  of  Equations  IV-9  and  IV-10  are  proportional  to  each 
other.  This  observation  leads  to  a  simple  ratio  for  the  time-parameters 
B^  and  B2  as: 

Bg/Bj  »  (IV-11) 

Consequently,  the  definition  for  the  Interface  parameter  x  can  be 
written  as: 
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2 

A  3  V(k2e2c2/k1plc1  tanh^'ff1^75j((|)  -  1  jJ/tanhFj  (IV-12) 

2 

where  ^  *  8^a  . 

Thus,  only  one  of  the  two  Equations  IV-9  and  IV-10  needs  to  be 
solved.  The  boundary  condition  of  ^  is,  as  usual, 

"Bi  -  e-  -*■  0 

By  defining  a  convenient  quantity  H  as: 

H  3  y-y-y  ItanHF'j ]/¥-|  (IV-13) 

Equation  IV-9  is  reduced  to: 

—  3  4HB?  ( IV- 14) 

d0  1 

Its  solution  with  boundary  condition  IV-13  gives  both  and  ^ 
variations  with  time.  Appendix  A  details  the  solution  method  of 
Equation  IV-14. 

Shown  in  Figure  IV-3  is  the  distribution  curve  for  c^/oj  *  3, 
kg/k-i  ■  10  and  b/a  *  2.  For  the  same  parametric  values,  the  -  IT 
curve  for  a  two-dimensional  equivalent  case  is  also  indicated  for 
comparison. 
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IV. 4  THE  AXIAL  DIFFUSION  EQUATIONS 


Re-expressinq  Equations  IV-7  and  IV-8  In  terms  of  T]  and  T2> 
one  has  the  following  equations: 


3T1  32T1  4H6^ 

=  “=r  +  n  -  H(1  +  R)]  [T2  " 

aT2  32T2  (k^kg)  4HI2 

^  7T  "  W  '  t(b/a)2  -”l  ]  '  ^  +  R»  l'2  "  h  J  UV‘,b; 

where  R,  the  thermal  capacity  ratio,  is  defined  by 

Pi  Ci 

R  3  - L-L— 5 -  ( IV- 17) 

P2c2[(b/ar  -  1] 


(IV-15) 


r  ^  ^  « 


/*(<  1  >-  \ 


In  parallel  with  the  developments  for  the  planar  case,  let 
4HB? 

L  ■  V~~h(V  +~rT  »v-'8) 

and  the  axial  diffusion  equations  can -be  represented  by  the  following, 
which  correspond  to  those  in  Section  III: 
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(IV-19) 


37,  32T, 

^V*l(VV 


(IV-20) 


The  inter-region  conductance  L  has  been  pre-determi ned  through  the 
l^-variation  obtained  from  Equation  IV-14.  See  Appendix  A  for  details. 

Because  of  the  close  parallel  developments  with  the  case  of  heat 
flow  along  laminated  composites,  the  solutions  and  the  methods  of 
obtaining  the  solutions  discussed  in  Section  III  are  directly 
applicable  to  the  above  equations  for  heat  flow  in  fibrous  composites. 
In  fact.  Equations  IV-19  and  IV-20  are  identical  to  Equations  III-19 
and  III-20,  with  the  only  differences  In  the  definitions  of  L  and  R 
(See  Equations  111-17  and  IV-17  for  R  and  Equations  III  18  and 
IV-18  for  L).  Numerical  results  obtained  for  composites  with 
uni -directional  fibers  are  discussed  in  Section  V. 
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V.  PERFORMANCE  EVALUATION  OF  THE  METHODS 


In  order  to  bring  out  the  essential  features  of  the  methods 
developed  in  this  report  as  well  as  significant  thermal  characteristics 
in  composite  bodies,  a  number  of  calculations  were  performed.  For 
simplicity  of  presentation  and  discussions,  only  composites  with 
uni -directional  cylindrical  fibers  were  analyzed  as  conclusions  and 
observations  deduced  for  this  configuration  are  equally  applicable  to 
laminated  composites. 

Principally,  two  sets  of  physical  parameters  were  employed  in 
the  demonstrative  calculations  presented  in  this  section.  They  were 

Specification  (A)  - 

k-|/k2  *  5,  ai/a2  s  3.33,  b/a  =  2,  R  =  0.5,  *  0.75 

Specification  (B)  - 

yk2  »  0.2,  ai/a2  -  0.2,  b/a  -  1.3,  R  »  1.45,  L.  -  10.5 

In  Specification  (A),  the  fiber  (material  1)  is  taken  to  be  more 
conductive  than  the  matrix,  indicative  of  applications  to  carbon/epoxy 
composites.  The  diffusivity  ratio  of  a 2 /o^  ■  0.3  is  to  allow  for  a 
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temperature  cross-over  in  the  Interior  of  the  composites  in  the  case 
of  constant  heat  flux.  In  Specification  (B),  reflective  of  metal 
matrix  applications,  the  fiber  conductivity  is  less  than  that  of 
the  matrix  material.  In  both  sets  of  the  specifications,  the  thermal 
capacity  ratio,  R,  is  In  the  order  of  unity,  l.e.,  the  two  media  have 
comparable  thermal  intertias.  The  asymptotic  values  of  the  transverse 
conductance  L  are  noted  to  be  L^,  a  0.75  and  10.5  respectively  for 
each  material.  It  turns  out  through  the  following  analysis  that  U> 
is  a  significant  quantity  in  composite  heat  transfer  analysis. 

Analyses  were  chiefly  focused  on  the  case  of  constant  surface 
heat  flux,  as  it  affords  a  more  severe  test  on  the  ability  of  the 
analytical  means,  and  more  importantly  it  represents  a  practical 
need  in  thermal  protection  technology. 


V.l  THE  ANALYTICAL  METHODS 

For  each  specification  of  the  physical  parameters,  several  methods 
of  analysis  were  employed  in  order  to  compare  their  results,  and 
thereby  to  assess  their  relative  merits  from  the  viewpoints  of 
accuracy  and  time-consumption.  The  methods  of  analysis  range  from 
the  most  elementary  approach  to  full  exact  solutions  of  the 
complete  diffusion  equations,  and  consist  of  the  following: 

(1)  Separate  Independent  solutions  which  assume  that  the 
fibers  and  the  matrix  material  have  no  transverse 
conduction.  In  this  situation,  each  region  behaves 
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Independently  of  the  other,  with  the  temperature  responses 
the  same  as  that  for  a  semi-infinite  domain  (x  >  0)  and  given 
by  the  following  expressions: 


erfc(x/2/ot.j8) 

(V-l) 

erfc(x/2*'a£F) 

(V-2) 

for  the  constant  surface- temperature  boundary  condition,  and 

T1  *  2(Q/kj  )»6p!F  1erfc(x/2v'Spr)  (V-3) 

T2  *  2(Q/k2)vEpF  ierfc(x/2/ct2&)  (V-4) 

for  the  constant  heat  flux  boundary  condition. 

(ii)  Double  heat-balance  integral -method,  which,  as  discussed 
in  Section  III,  treats  the  Inter-region  transverse  heat 
conduction  on  an  integrated  basis.  The  axial  variation  of 
the  temperature  Is  characterized  by  a  single  parameter  Yj 
or  y2  In  each  material  region.  Variations  of  the 
characteristic  time-dependent  parameters  are  determined 
by  integrating  (with  respect  to  8)  Equations  III-27  and 
III-28. 

(Ill)  Semi -finite  difference  solution  which  Is,  in  essence,  a 
finite-difference  scheme  for  solving  the  axial  diffusion 


Equations  I I 1-1 9  and  I I 1-20.  Since  the  transverse 
temperature  variation  Is  expressed  through  the  time- 
dependent  parameter  6^,  the  use  of  a  finite-difference 
to  account  for  the  axial  temperatur  <tr1at1on  leads  to 
the  terminology:  semi -finite-difference  solution. 

(1v)  Exact  solution  or  finite-difference  solution  of  the 
complete  diffusion  Equations  III-8  and  III-9,  while  it 
Is  not  the  objective  of  the  present  investigation  to  analyze 
heat  conduction  in  composites  by  exact  solutions,  results 
thus  obtained  serve  as  a  reference  by  which  other  methods 
can  be  judged.  Because  of  time  and  storage  requirements, 
only  a  limited  number  of  computations  were  made,  as  will 
be  noted  later  in  this  presentation. 


V.2  TEMPERATURE  DISTRIBUTIONS  FOR  SPECIFICATION  (A) 

For  the  parametric  values  in  Specification  (A),  the  transverse 
time-dependent  parameter  6^  Is  first  determined  by  the  numerical  method 
outlined  In  Appendix  A.  Figure  V.l  shows  how  8^  varies  with  ef,  the 
non-dimensional  time.  Also  Included  In  Figure  V.l  Is  the  transverse 
conductance  l,  which  reaches  an  asymptotical  value  of  l  *  0.75  at 
large  values  of  e.  In  practice,  however,  the  transverse  conductance 
remains  virtually  unchanged  at  L  »  0.75  for  "5T  ■  1  or  larger.  Physically, 
It  means  that  the  transverse  temperature  profiles  in  the  fiber  and 
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matrix  regions  have  reached  their  equilibrium  shapes,  though  their 
magnitudes  may  still  be  changing  owing  to  axial  diffusion.  The  time 

for  such  an  event  is  characterized  by  0  s  1  or  the  physical  time 

2 

6  =  a  /Oj,  the  time  necessary  to  diffuse  over  a  distance  of  the  fiber 
or  tow  radius  of  a. 

As  the  boundary  condition  of  constant  surface  heat  flux  (equal  in 
both  regions)  represents  a  more  practical  situation  encountered,  the 
surface  temperature-rise  with  time,  as  heating  proceeds,  is  a 
significant  measure  of  the  material  capability  to  diffuse  the  heat 
flux  away  from  the  surface.  Figure  V.2  shows  the  matrix  surface  and 
fiber  surface  temperature  variations  according  to  different  methods 
of  analysis. 

First,  Method  (i),  which  assumes  that  each  region  responds 
independently  of  its  adjoint  neighbor,  produces  the  extremes  of  the 
temperature- time  curves.  For  the  physical  parameters  specified,  the 
matrix  surface  has  a  faster-rising  temperature  than  does  the  fiber 
surface.  Since  any  mutual  conduction  in  the  transverse  direction 
tends  to  bring  these  two  temperature  responses  closer  together,  the 
calculated  temperature  responses  by  Method  (i)  therefore  constitute 
an  outer-limit  envelope. 

Conversely,  the  use  of  full  diffusion  equations,  i.e..  Method 
(iv),  produces  the  exact  solutions  to  the  problem,  where  the  Inter¬ 
region  conduction  is  exactly  accounted  for,  constitutes  an  inner- 
limit  envelope.  Between  the  inner  limit  and  the  outer  envelope, 
there  exists  a  fairly  wide  latitude  in  which  the  results  of  the  other 
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two  methods  lie.  Results  based  on  the  semi -finite-difference  method, 
i.e.  Method  (iii),  are  located  quite  closely  to  the  exact-solution 
results,  whereas  results  based  on  the  double  heat-balance  integral  - 
method,  i.e.  Method  (iv),  are  closer  to  those  from  Method  (i)  than 
other  methods. 

Proximity  to  the  exact-solution  results  appears  to  lend 
credence  to  Method  (iii),  i.e.,  the  semi -finite-difference  solution 
of  Equations  III-19  and  III-20,  although  further  improvements  are 
desirable  and  indeed  possible  through  a  better  representation  of 
the  transverse  temperature  distributions  in  both  regions.  In¬ 
ability  of  Method  (ii)  to  include  transverse  conduction  effect  is 
attributable  to  its  one-parameter  (y^  or  Yg)  representation  of  the 
axial  temperature  profile,  and  to  the  fact  that  inter-region  conduction 
is  only  accounted  for  on  an  integrated  basis.  Along  the  axial 
direction,  cross-over  of  the  temperature  profiles,  as  will  be 
discussed  next,  indicates  reversal  of  transverse  heat  conduction 
from  one  direction  to  another.  Thus,  in  summing  up  the  local 
transverse  conduction  distribution,  mutual  cancellation  causes 
significant  errors  in  the  final  results  of  analysis. 

To  further  examine  the  relative  merits  of  the  various  methods 
of  analysis,  the  temperature  variations  along  the  axial  direction  are 
shown  in  Figure  V.3.  Only  the  results  from  Methods  (iii)  and  (iv) 
are  presented  and  the  closeness  between  the  results  from  Method  (iii) 
to  those  from  the  exact  method  gives  another  indication  of  the 
usefulness  of  the  semi -finite  difference  method.  Figure  V.3  also 
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reveals  as  to  why  the  double  heat-balance  integral -method  fails  to 
yield  reasonable  results  because  of  the  cross-over  of  the  temperature 
profiles.  For  larger  times,  W  =  1  and  W  s  2,  the  temperature  profiles 
from  Method  (iii)  are  compared  with  those  from  Method  (i),  i.e.,  the 
no-interaction  solutions  in  Figures  V.4  and  V.5.  A  significant 
observation  based  on  the  data  in  these  three  figures  is  the  location  of 
the  cross-over  point  at  (x/a)  t  1.0;  its  precise  position  varies  not 
appreciably  from  the  value  of  (x/a)  3  1  in  these  figures.  The  near 
constancy  of  the  cross-over  location  is  in  consonance  with  the  diffusion 
of  ¥  =  1  or  a  diffusion  distance  of  one  fiber  radius  a.  Hence  in 
problems  involving  constant  surface  heat  flux,  large  temperature 
disparity  in  two  different  media  in  thermal  contact  is  expected  to  be 
confined  to  a  surface  layer  of  no  more  than  a  few  fiber  or  tow  radii 
deep.  Temperature  differentiation  outside  this  layer  is  substantially 
reduced  by  mutual  transverse  conduction.  From  a  thermal  stress  viewpoint, 
it  is  in  this  surface  layer  where  high  stress  values  would  be  found. 

As  a  further  confirmation  of  the  semi -finite-difference  approach, 
the  case  of  a  step-rise  in  surface  temperature  is  also  investigated. 

Figure  V.6  contains  a  typical  record  of  the  calculated  distributions 
and  those  of  the  exact  solutions.  It  becomes  apparent  that  the  semi- 
finite-difference  method  gives  excellent  agreement  with  the  exact 
solutions.  Comparing  the  results  shown  in  Figure  V.3  with  those  in 
Figure  V.6,  it  should  be  noted  that  the  larger  difference  exists  in 
Figure  V.3  than  in  Figure  V.6.  For  the  former,  the  boundary  condition 
is  that  of  constant  heat  flux;  for  the  latter  it  is  that  of  a  constant 
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Figure  V.4.  Comparison  of  Axial  Temperature  Distributions  in  Fiber 
and  Matrix  0=1.  For  Constant  Heat  Flux  (k./k,  =  5, 


surface  temperature.  It  becomes  understandable  if  it  is  re-called 
that  the  solution  to  a  constant  surface-heat  boundary  condition  can  be 
synthesized  to  consist  of  an  infinite  number  of  small  surface 
temperature  rises.  In  this  fashion,  the  minute  errors  as  shown  in 
Figure  V.6  are  accumulated  to  the  error  exhibited  in  Figure  V.3. 

V.3  TEMPERATURE  DISTRIBUTIONS  FOR  SPECIFICATION  (B) 

The  physical  parameters  in  Specification  (B)  were  selected  with 
metal  matrixes  in  mind  in  which  the  fibers  have  lower  conductivities 
than  the  matrix  materials.  A  thermal  diffusivity  ratio  of  0.2  was 
assumed  so  that  cross-over  of  the  axial  temperature  profiles  for  the 
two  media  would  occur  for  constant  surface-heating  boundary  condition. 
As  a  large  number  of  composites  in  aero-space  applications  have  a 
volume  ratio  of  0.6,  a  diameter  ratio  of  1.3  was  considered  for  (b/a) 
which  yields  a  volume  ratio  of  0.69.  The  preceding  parameters  result 
in  a  heat  capacity  ratio  of  R  a  1.45  which  assures  that  the  axial 
temperature  distributions  are  substantially  modified  from  their 
independent  (no  transverse  conduction)  solutions. 

Shown  in  Figure  V.7  are  the  variations  of  the  transverse  time- 
parameter  Uj  and  the  transverse  conductance  L.  The  latter  reaches  an 
asymptotic  value  of  *  10.4,  a  value  substantially  higher  than  that 
for  Specification  (A),  but  in  a  much  shorter  time.  Consequently,  the 
surface  temperature  variations  (with  time)  for  the  two  materials  show 


a  much  narrower  spread  as  Indicated  by  the  data  In  Figure  V.8  than  those 
based  on  Specification  (A),  wherein  the  asymptotic  transverse  conduc¬ 
tance  is  only  0.75.  Hence,  It  can  be  concluded  that  a  governing 
parameter  Is  the  transverse  conductance  L.  A  larger  value  of  L 
signifies  effective  transverse  conductance  which  brings  about  equalization 
of  the  temperature  in  the  two  regions  in  a  much  shorter  distance  from 
the  end  where  heating  starts.  Consistent  with  the  arguments  set  forth 
in  the  preceding  are  the  temperature  responses  shown  in  Figure  V.9 
where  the  cross-over  of  the  temperatures  in  the  two  regions  takes  place 
near  (x/a)  =  0.5  to  0.6  for  0"  =  0.4  and  1.0  respectively.  Beyond 
the  cross-over  point,  the  axial  temperature  distributions  are  nearly 
parallel  to  each  other  with  a  much  smaller  differential  between 
those  shown  in  Figures  V.3  and  V.4,  for  which  a  distinguishing  feature 
is  a  lower  transverse  conductance  of  L  =  0.75. 

oo 

Examination  of  the  formula  which  defines  the  asymptotic  value  of 
the  transverse  conductance,  L^,  indicates  that  the  make-up  parameter 
Is  the  fiber/matrix  conductivity  ratio  modified  by  the  volume  ratio 
of  the  two  materials.  As  .a  further  indication  of  the  significant 
role  of  this  ratio,  the  transverse  temperature  distributions  across 
the  fiber  and  the  matrix  regions  are  shown  In  Fgiure  V.10  based  on 
the  data  for  Specification  (B).  The  fiber  has  much  less  conduction 
than  the  matrix  (k1/k2  =  0.2),  resulting  in  a  wider  temperature 
variance  in  the  fiber  region  than  in  the  matrix  region.  At  (Va)  *  0, 
i.e.,  on  the  heating  surface,  and  for  F  *  0.4,  the  fiber  cross-section 
has  a  very  substantial  temperature  variation  from  the  center  to  its 
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Figure  V.8.  Surface  Temperature  Responses  to  Impulsive  Heat  Flux, 
Fiber-Composite  (k^/k^  *  0.2,  dj/cx^  *  0.2,  b/a  *  1.3). 
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Figure  V.9.  Axial  Temperature  D1 
0  3  0.4  and  ¥  3  1  fo 
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Heating  Surface 


Figure  V. 10.  Transverse  Temperature  Distributions  at  Various  Axial 
Positions,  8  -  0.4  for  Impulsive  Surface  Heat  Flux, 
(Icj/kg  8  0.2,  ot.|/a2  »  0.2,  b/a  ■  1.3). 


edge.  High  Transverse  conductance  rapidly  equalizes  the  two  regions 
at  (x/a)  >  0.4,  even  though  relative  variations  in  these  two  materials 
still  exist  In  the  Interior  of  the  composite. 


V.4  LARGE  THERMAL  CAPACITY  RATIOS 

As  an  additional  exploration  of  relevant  parameters  In  transient  flow 
in  composites,  several  computer  runs  were  made  to  decipher  the  effect  of 
the  thermal  capacity  ratio  R  on  the  temperature  distribution.  A  typical 
record  is  illustrated  in  Figure  V.11  for  k-|/k2  3  10.  A  thermal  diffuslvity 
ratio  (oj/ct^)  of  unity  Is  taken.  These  parameters  and  others  resulted  in 
a  thermal  capacity  ratio  of  R  3  14.5,  which  expresses  the  condition  that 
the  fiber  has  a  much  larger  thermal  Inertia  characterized  by  the  product 
of  (pc)  and  the  fiber-volume  than  that  of  the  matrix  region.  The  fact 
that  ctj/cig  ■  1  precludes  a  temperature  cross-over,  as  Is  evident  in 
preceding  combinations.  The  temperature  distribution  curves  for 
¥  3  0.4  and  15-  3  1  show  a  practical  merge  of  the  temperatures  In  the 
two  different  media  at  (x/a)  3  1,  which  would  become  the  cross-over 
point  If  <  ctg. 

As  a  reference  for  discussion,  the  fiber  temperature  distribution 
for  S’  3  0.4  by  Method  (1)  Is  also  Included  In  Figure  V.ll.  This 
reference  distribution  Is  of  course  based  on  the  condition  that  the 
fiber-region  temperature  Is  Independent  of  that  In  the  matrix  region. 

The  corresponding  temperature  distribution  for  the  matrix  region  Is  not 
present  In  Figure  V.ll,  for  It  would  be  situated  above  the  scope  of  the 
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Figure  V.ll.  Temperature  £1str1but1ons  of  Fiber  and  Matrix  at 
9  *  0.4  and  W  =  1  for  Impulsive  Heat  Flux. 
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ordinate.  The  influence  of  the  thermal  capacity  R  is  accordingly 
apparent:  because  of  a  mutual  thermal  interaction  by  virtue  of  transverse 
conduct! on, the  fiber  temperature  is  raised  and  the  matrix  temperature 
is  lowered  from  their  respective  reference  distribution  curves  for 
zero  transverse  condition.  The  deviation  from  the  reference  curve 
is  of  course  more  for  the  matrix  region  than  for  the  fiber  region, 
but  not  in  exact  proportion  to  the  value  R.  Axial  diffusion  modifies, 
to  some  extent,  the  transverse  shifts  of  the  temperature  distribution 
curves. 
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VI.  CONCLUSIONS 


The  analyses  presented  in  this  report  consider  the  transient 
heat  flow  along  uni -directional  fibers  of  composites  or  in  the  stacking 
plane  of  laminated  composites.  Because  of  the  preponderance  of 
physical  and  geometrical  parameters  involved,  exact  solutions  by 
analytical  means  are  usually  not  feasible,  except  for  special  cases 
of  very  much  simplified  geometries.  Furthermore,  sorting  through 
a  large  number  of  exact  numerical  solutions  in  order  to  extract 
important  parametric  groups  is  prohibitively  difficult. 

Here  in  this  report,  a  method  has  been  developed  which 
identifies  the  parametric  combinations  and  is  relatively  amenable 
to  numerical  solutions.  The  method  has  two  components:  (1)  the 
transverse  (to  fibers  or  laminations)  temperature  profiles  are 
treated  by  means  of  two  transverse  time-dependent  parameters  which 
are  solutions  to  heat-balance  differential  equations,  one  for  each 
region.  (2)  The  resulting  equations  for  the  axial  temperature 
variations  in  the  two  materials  are  much  simplified,  and  their 
numerical  solutions  are  greatly  simpler  to  implement  than  the  full 
diffusion  equations.  The  axial  diffusion  equations  contain  the 
governing  parametrical  groups  which  are  identified  to  play 
Important  roles  in  the  temperature  responses  to  a  heat  flux  at  an 
exposed  surface.  They  are  (1)  a  transverse  conductance  parameter  and 
(11)  the  thermal  capacity  ratio  of  the  two  different  materials 
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comprising  a  composite.  The  first  governs,  for  the  case  of  constant 
surface  heating,  the  extent  of  a  depth  near  the  heating  surface  wherein 
a  large  temperature  differential  exists  between  the  two  materials.  The 
second  parameter  influences  the  relative  temperature  modifications  from 
their  respective  reference  values,  the  latter  being  those  based  on  the 
idealized  one-domain  no-interaction  solution  for  each  region  separately. 
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APPENDIX  A.  NUMERICAL  SOLUTION  OF  THE  TRANSVERSE 
HEAT-BALANCE-INTEGRAL  EQUATION. 


In  the  transverse  direction,  thermal  diffusion  time  -  time 
required  to  even  out  the  transverse  temperature  variations  in  the  two 
different  media  in  contact  -  is  much  smaller  than  thermal  diffusion 
time  in  the  axial  direction.  This  is  due  primarily  to  a  difference  in 
the  dimensions  along  these  directions  in  a  composite  material  with 
thin  laminations  or  fibers  of  micro-dimensions  compared  to  the  length. 
The  transverse  heat  flow  can  therefore  be  treated  independently  of  the 
overall  heat  conduction  process  and  this  is  accomplished  by  a  heat- 
balance  integral  method  as  outlined  in  Sections  II  and  III  for  two- 
dimensional  heat  flow  (laminated  composites)  and  uni -directional 
f  i  brous-composi  tes . 

The  resulting  Equations  III-8  and  IV-9  are  grouped  together  as: 


d_ 

d0 


^  tanhBj 


1  +  X  T 


31  J 


1  nx 

"  T+T 


61tanhB1 


(A-l) 


X 


*C2P2C2  tanh^l 
klplCl 


tanhB^ 


(A— 2 ) 


where  the  factors  n  and  r  assume  the  following  values: 

n  *  4,  r  »  2  for  axis-symmetrical  heat  flow, 

n  *  1,  r  *  1  for  two-dimensional  heat  flow. 

The  boundary  condition  on  ^  Is  ^  as  e  -*■  0. 

As  o’  •*  0,  •*> «,  the  parameter  A.  assumes  the  following  value: 


2P2  2J 


MP1*T 


Equation  A-l  is  reduced,  for  small  vaues  of  0,  to: 


(A-3) 


(A— 4) 


for  which  the  solution  is  clearly: 


(A- 5) 


For  large  values  of  0,  the  asymptotic  solution  A-5  breaks  down 
and  the  full  Equation  A-1  must  be  used  to  obtain  the  solution.  To  do 
this,  the  equation  is  cast  in  the  form: 
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g.  jLUL  d  [r^-r  -"hBl 

nXg^tanhe^  1  +  X  81 


Or,  in  terms  of  the  parameter  H  defined  in  Section  III,  the  above 
equation  is  put  as: 


d?  *  Vd(H) 


(A-6) 


The  parameters  V  and  H  are  reproduced  here  for  clarity: 


V  = 


1  +  X 


nXB1  tanhB-j 


(A- 7) 


^  tanhBj 


1 


(A-8) 


By  assigning  a  series  of  values  to  ftj,  starting  from  a  large 
value  and  down  to,  say,  7^  *  0.1,  a  curve  In  the  H-V  plane  is  obtained. 

From  Equation  A-6,  it  is  clear  that  9"  represents  the  area 
underneath  the  H-V  curve  as  sketched  in  Figure  A.l. 

As  -*■  0,  H  and  X  approach  their  limits  as  follows: 


X  ■  A*  •  (pgCg/p^CjJKb/a)1*  -  1]  «  1/R 


(A-9) 


H  *  Hob  »  Xco/O  +  Xoo) 


89 


Hence  for  very  large  values  of  6,  the  numerical  procedure  Is 
Inherently  Inaccurate  as  H  -*•  H,,  and  V  Increases  without  limit.  For 
this  range,  l.e.,  ^  -*■  0,  the  parameters  H  and  V  are  expressed  In 
their  asymptotic  forms  as: 


(A-10) 


(A-1 1 ) 


where  s  *  [1  +  (k1/k2)((b/a)r  -  1)]/3[1  +  1/A.l. 

Substitution  of  these  asymptotic  forms  into  Equation  A-6  yields: 


*  g1  exp  [-  (n/s)e  ]  (A-12) 

The  constant  of  Integration  g^  can  be  determined  by  fitting  the  above 
equation  to  the  last  point  of  9^  ~  8^  curve  obtained  from  the  numerical 


procedure. 


APPENDIX  B 

FINITE-DIFFERENCE  SOLUTION  OF  THE  AXIAL 
HEAT  DIFFUSION  EQUATIONS 
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APPENDIX  B.  FINITE-DIFFERENCE  SOLUTION  OF  THE  AXIAL 
HEAT  DIFFUSION  EQUATIONS. 


With  temperature  variations  In  the  transverse  direction 
approximated  by  two  form  functions  with  two  transverse  time-parameters 
81  and  &2*  the  reduced  heat  conduction  equations  in  the  axial  direction 
are  those  of  Equations  III-19  and  I I 1-20  for  the  two-dimensional  case 
(laminated  composites)  and  Equations  IV-19  and  IV-20  for  the  axis- 
symmetrical  case  (fibrous  composites).  These  two  sets  of  equations 
can  be  jointly  represented  by  the  following: 


o2T, 

~  Tl} 


«2  o2T2 
a1  1? 


-  RL(T2  -  T}) 


(B-l) 


(B— 2 ) 


The  Interaction  parameter  L  denotes  the  Inter-region  heat  conduction 
effects  based  on  their  respective  average  temperatures  over  each 
cross-section  and  Is  defined  by: 


The  thermal  capacity  ratio  R  Is  defined  by: 


R  -  (P1c1/p2c2)/[(b/a)j  -  1]  (B-4) 

Equation  B-3  for  L  and  Equation  B-4  for  R  cover  the  two  cases  through 
the  coefficients  n  and  j  In  these  two  equations.  They  have  the 
following  numerical  values: 

(1)  For  two-dimensional  heat  flow  (laminated  composites) 

n  a  1 ,  j  =  1 

(il)  For  axis-symmetrical  heat  flow  (fiber-composites) 

n  ■  4,  j  *  2 


ASYMPTOTIC  ANALYSIS  FOR  SMALL  TIMES.  9^0 

It  is  of  Interest  to  note  that  at  very  small  times,  the  transverse 
conduction  parameter  exhibits  a  singularity  for  an  Impulsive  type  of 
thermal  loading,  either  an  Impulsive  temperature  rise  or  an  impulsive 
heat  flux  at  x  »  0.  The  analysis  In  Appendix  A  shows  that  ^  -  l/VT, 
which  Is  conslstant  with  more  rigorous  theoretical  results.  As  the 
intermediate  parameters,  H  and  X,  are  defined  by: 
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H  -  [tanhS,  1/5,  (B-5) 


the  small-time  variation  of  L  can  be  determined  as  follows: 

$  -*■  0,  S,  ■  1/V2nf,  X  XQ 

where  XQ  *  v^l^P^cpTTI^P^c^T  *  (kg/k-,  )vS^7o£*  Hence,  the  intermediate 
parameter  H  becomes: 

H  -  tXg/O  +  XQ)lV2ne 

and  the  inter-region  conduction  parameter  L  is  of  the  form: 


(B— 7) 


ASYMPTOTIC  ANALYSIS  FOR  LARGE  TIMES,  9  -  «° 

As  the  inter-region  conduction  proceeds,  the  interaction  parameter 
L  reaches  a  " steady- state"  value,  while  the  temperature  profile  in  each 
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region  establishes  Itself  In  a  steady-state  pattern  In  the  transverse 
direction.  Under  this  condition,  6"  -►  «,  fTj  0,  and  the  Intermediate 
parameters  H  and  X  are  analyzed  below: 

9^  °°»  8^  3  g1exp[-(n/s)0]  (B-8) 

which  Is  obtained  In  Appendix  A  [Equation  A-12],  where: 

s  -  (1  +  (k1/k2)[(b/a)j  -  1  ]}/3[1  +  R]  (B-9) 

Other  relevant  parameters  can  be  expressed  as: 


(B-9) 


*-{' -r[^  (<#>J  - 1  )Z ' 1  ]  *  -  •  j 

(B— 1 0) 

HJ\  -  s  6^  +  .  .  .  ] 

(B— 1 1 ) 

LJ\  -  si?  +  .  .  •] 

(B-12) 

(1/R)  -  (p^/p^^Kb/ajj  -  1] 

(B-13) 

1/(1  +  R) 

(B-14) 

3n/{1  +  (k1/k2)[(b/a)J  -  1]} 

(B-15) 

1 


The  last  quantity  Is  of  particular  Interest,  for  It  represents  the 
"steady- state"  Inter-region  conductance,  which  was  analyzed  by  Maewal, 
Gurtman  and  Hegemeir  (J.  of  Heat  Transfer,  vol.  100,  p.128,  Feb.  1978). 
In  their  analysis,  the  variation  of  L  at  small  iT-values  was  not  taken 
into  account  for  a  first-order  estimate.  Figure  B-l  illustrates  the 
variation  of  L  as  a  function  of  time  0"  for  a  combination  of  k?/k1  =  10, 
a2/a-|  *  3  and  b/a  =  2  and  for  the  axis-symmetrical  case  (fibrous 
compos i tes ) . 


FINITE-DIFFERENCE  EQUATIONS 

Let  the  iC-axis  be  divided  into  a  number  of  equi-distant  nodes; 
starting  from  x  =  0,  these  nodes  are  numbered  by  an  index  j  =  1,  2, 
3,...  .  The  distance  between  two  adjacent  nodes  is  AX'  and  the  time 
advance  is  represented  by  A6\  A  three-point  explicit  marching  (in 
time)  scheme  is  used, with  the  following  diagram  illustrating  a 
three-point  cluster  of  nodes  (j  -  1),  j,  and  (j  +  1): 


Ti.  J-l  T1.J  Tl,  J+l 

— • - • . -  — » 

j-l  j  j+l 

•*  x" 
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Equation  B-1  is  reduced  to: 


ti,j  *  Tu  .  i 


a?  (ajO 


?  [Ti,  j+i  +  Ti,  j-i  '  2TU] 


+  LtT2.J  -  Tl,jl 


(B-1 6) 


where  the  unprimed  quantities  are  current  values  and  the  primed 

_ I 

quantity  T,  .  is  the  future  temperature.  Casting  B-1 6  as: 

*  f  J 


T,  ..  =•  [T 


+  T,  ,,)  +  T.  ,[1  -  2 


i f j  '^jfr  l,  j+i  i.  j-iJ  i,j 


A0 


(A*2) 


LAB] 


+  LA0[T2,j’ 


( B-1 7 ) 


A  similar  procedure  on  Equation  B-2  yields: 


j  x  w2  AB  rT  .  T  , 

T2,j  *2,  j+1  T2,  j-lJ 

+  T.  4  [  1  -  2-^f  -  RLAB]  +  RLA8[T,  .]  (B-18) 

2’j  (AX)Z\°1/  1J 

Stability  of  the  explicit  marching  scheme  Is  assured  by  choslng  AB 
such  that  the  coefficients  of  T1  j  and  T2  j  on  the  right-hand  sides  of 
B-1 7  and  B-18  remain  positive. 
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APPENDIX  C 

EXACT  SOLUTIONS  FOR  LAMINATED  COMPOSITES 
OF  FINITE  LENGTH 
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APPENDIX  C.  EXACT  SOLUTIONS  FOR  LAMINATED  COMPOSITES 
OF  FINITE  LENGTH 


C-I.  DEVELOPMENT  OF  SOLUTIONS 

Consider  the  configuration  depicted  in  Figure  C.l  which  shows  a 
repeating  section  of  a  laminated  composite  stacked  by  alternate 
layers  of  two  different  materials  1,  and  2.  The  extent  in  the  x- 
direction  is  finite. 

At  Xs  L,  where  heating  occurs,  two  boundary  conditions  are 
considered:  one  is  for  an  impulsive  temperature  rise  from  an  initial 
temperature  distribution  -  taken  to  be  zero;  the  other  is  for  an 
impulsive  heat  input  -  equal  in  both  regions.  These  two  boundary 
conditions  are  the  basic  ingredients  with  which  solutions  of  more 
practical  interest  can  be  built  by  a  superposition  principle. 
Designating  these  two  boundary  conditions  as  the  T-case  and  Q-case 
respectively,  these  can  be  expressed  by: 

T-case  Tj  »  1 ,  Tg  *  1 ,  at  X  *  L  (C-l ) 

Q-case  kix*Vax"  k^tyax  -  Q,  at  X  «  L  (C-2) 
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Figure  c.?. 
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At  the  other  end  x  «  0,  two  types  of  boundary  conditions  are  considered; 
they  are: 

(i)  T-j  *  Tg  *  0  at  x  »  0  (Fixed  Temperature)  (C-3a) 

(ii)  3^/ax  *  STg/Sx  *  0  at  x  *  0  (Insulated)  (C-3b) 

The  initial  condition  is  of  course  a  uniform  (taken  to  be  zero) 
temperature  distribution.  At  different  y-positlons,  the  following 
conditions  apply: 

3Tj/3y  =  3T2/3y  =0  at  y  *  0,  and  y  =  b  (C-4) 

T,  «  T2  and  k^s^/sy  *  k2y3T2/3y  at  y  -  a  (C-5) 

All  told  therefore,  four  cases  are  treated  In  this  report;  they 

are: 


Case  (1):  T^y.L.e)  ■  1, 

T2(y,L,8)  •  1 

(C-6a) 

aT^y.o.eJ/ax  - 

0,  3T2(y,0,e)/3x  ■  0 

(C-6b) 

Case  (2):  Tj(y,L,e)  ■  1, 

T2(y.U8)  •  i 

(C-7a) 

T^y.o.e)  •  o. 

T2(y,0,e)  ■  0 

(C-7b) 
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Case  (3):  ^lxC3T1 (y,8)/3x]  ■  Q,  k2x[3T2(y,e)/3X]  »  Q  (C-8a) 

aT^y.o.eJ/ax  ■  o,  3T2(y,o,e)/ax  «  o  (c-8b) 

Case  (4):  k^taT^y.L.eJ/ax]  *  Q,  k2x[3T2(y,L,0)/3x]  *  Q  (C-9a) 

^(y.0,6)  ■  0,  T2(y,0,9)  -  0  (C-9b) 


All  these  cases  are  subjected  to  the  interface  boundary  conditions 
specified  by  Equations  C-4  and  C-5  plus  the  initial  conditions  of: 

T^y.x.O)  -  T2(y,x,0)  *  0  (C-10) 

Before  procedlng  with  the  Individual  cases,  the  governing 
differential  equations  are  first  non-dimensional  1 zed.  Let  the 
characteristic  length  be  L  and  the  characteristic  time  be  (L  /<*|X). 
the  following  definitions  and  variable  transformations  are  obtained: 

°ix  *  kix/(pc)r  ^  *  k2x/(pc)2*  y  *  y/L.  *  *  x/0 

>  (C-ll) 

A1  *  (kl/klx>»  X2  •  (k2y/k2x)  J 

9  •  (o^e/l2),  a  ■  a/L,  E  -  b/L 

The  non-dimensional  equations  governing  the  temperature  distributions 
Tj ,  and  T2  for  these  two  regions  are: 
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(a^/ax2)  +  x1  (a^/ay2)  -  (aiyae)  (c-12) 

(a2T2/ax2)  +  x202T2/ay2)  =  (a^ctgJCaTg/ae)  (c-13) 


CASE  (1).  Constant  Temperature  Rise  at  One  End  and 
Insulated  Surface  at  the  Other 


The  complete  solutions  for  T1  and  T2  can  be  immediately  written  as: 


T 


1 


+  i  L  Ea,  e  8^j9[cosCiTmx/2):Y1(y) 
J-1,2..  m-1,3,5..  ^  1 


(C-14) 


T2 


2  _ 

'kz...Lz. 5.>  •  "* 


(C- 15) 


where  the  appearance  of  cos(imrx/2)  with  m  as  any  odd  integer  Is 
necessitated  by  Equations  C-6a  and  C-6b  at  7  *  1  and  3T  ■  0  respectively. 

Direct  substitutions  of  Equations  C-14  and  C-15  Into  Equations 
C-12  and  C-13  give  the  following  which  govern  and  72. 

+  [B2  -  (nrn/2)2]Y1  •  0  (C-16) 
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*2^2  +  “  ( nnr/2)ii]T2  ■  0 


(C-17) 


tquatlons  C-l 6  and  C-17  are  subjected  to  the  boundary  conditions 
(dYj/dy)  *  0  at  y“=  0,  (dYg/dy)  ■  0  at  y  =  F.  Furthermore,  because 
of  the  forms  of  Equations  C-14  and  C-15  and  the  equal  heat  flux 
requirement,  two  additional  boundary  conditions  appear: 


Y1 (a)  *  Y2(a) 
kly*l (a)  3  k2yT2(a) 


(C-l 8) 

(C-l 9) 


The  fact  that  there  are  two  homogeneous  differential  equations  to 
comply  with  four  homogeneous  boundary  conditions  results  in  an 
eigen-value  problem  for  8.  Specifically,  one  first  solves  for 
Y1  and  Y2  from  Equations  (C-l 6)  and  (C-17)  giving  the  following. 


l.mj 


r2,mj 


_  cosy ^[8*.  -  .(nm/2)2l/X1 

■ 

■  (imr/2)2J/A1 

cos(y  -  b)yj 

l(al/a2)e«d  ' 

cos  (a  -  F)^ 

ll(a l/a2)em  *  (nnr/2)2]/A2 

(C-20a) 


(C-20b) 


The  forms  of  Equations  C-20a  and  C-20b  have  been  adjusted  to  comply 
with  all  boundary  conditions  at  various  y's,  except  at  the  Interface 
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y  ■  a  which  supplies  the  eigen-value  equation  for  6^.  Performing 
the  operation  Indicated  In  Equation  C-5  gives  rise  to: 


(C-21 ) 

Equation  C-21  determines,  for  each  odd  integer  of  m,  a  series  of  B  .- 

mj 

values,  with  the  running  indexes  j  a  1,  2,  3,  ...  to  denote  successive 

roots.  The  eigen-value  Equation  C-21  is  discussed  in  Section  C-IV. 

For  each  eigen-value  B_,-  thus  obtained,  an  eigen  function  Y  .  is 

mj  mj 

therefore  defined  for  the  entire  domain  from  y  =  0  to  y  «  b.  The 
function  7^  is,  however,  defined  by  two  sub-region  functions  71  ^ , 
for  y  »  0  to  a,  given  by  Equation  C-13,  and  ^  mj*  for  7  *  to  b, 
given  by  Equation  C-20.  The  eigen-functions  7^  and  7m^(i  f  j)  are 
orthogonal  to  each  other  In  the  entire  domain  with  respect  to  a 
weight  function  w^  as  follows: 

y  ■  0  to  a  w^  ■  1 

(C-22) 

y  ■  a  to  F  wf  ■  (pc)2/(pc)1 
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The  development  of  the  preceding  steps  is  contained  in  Section  C-ll  and 
parallels  with  the  original  derivation  by  Title  and  Johnson  (ASME 
paper  No.  65-WA/HT-52). 

To  determine  the  coefficients  E^,  Equations  C-14  and  C-15 
are  combined  to  read: 


2  — 

00  00  .0 

T  *  1  +  £  L  E  .e  mj  cos(rmrx/2)Y  .  (C-23) 

j»l ,2, . .  m-1,3,5,..  ^  mj 

At  ¥  *  0,  the  Initial  temperature  distribution  is  zero  everywhere. 
Setting  T  to  zero  on  the  left  and  ?  to  zero  on  the  right.  Equation 
C-23  is  multiplied  by  [cos{rnnT/2)7n.w^]  and  is  integrated  from  y  =  0 
to  y  *  ¥  and  from  x  ■  0  to  x  3  1 .  The  resulting  expression  for  E^ 
turns  out  to  be  (see  Section  C-III  for  details): 

m-1 

E^  -  -(8/imr)(-l)_r[ft/(fa  +  ft)]  (C-24) 
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CASE  (2).  Constant  Temperature  Rise  at  One  End  and 

Temperature  Maintained  at  the  OtTiirTna — - 


The  solutions  in  regions  (1)  and  (2)  can  be  shown  as  follows: 


eLe 


T1  a  x  +  £  T.  F  .e  sinfmrxJY, 

'  nal ,2,3, . .  J-1,2,..  nj  1 


(C-25) 


V*+£  £  F  e^sinfmrx^ 

£  n-1 ,2,3, . .  J-1,2,..  nj  2 


(C-26) 


7  -7  -  C°Sy  V  [gn.l  - 

1  1 » n j  _  j  o  "9" 

COSaVie nj  '  ("»)  ]/*! 


(C-27) 


for  y  *  0  to  a,  and 


B  J  m  cos(y  -  b)\/ l(o1/tt2)B^1  -  (mr)Z]/X 

-  (nw)2]/X; 


r2  T2,nj 


COV.d 


(C-28) 


for  y  »  a  to  F. 
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Equality  of  the  Interface  temperatures  at  y  *  a  is  assured  by  the 
forms  of  Equations  C-27  and  C-28  and  that  the  same  coefficient  Fpj 
appears  in  Equations  C-25  and  C-26.  The  condition  of  equal  heat  flux 
on  both  sides  of  the  interface  yields  an  eigen-equation  for  Sn j : 


Equation  C-29  serves  to  determine,  for  each  integer  n,  a  series  of 
6pj-values  for  j  a  1,  2,  3,  ...  to  denote  the  successive  roots. 

The  coefficients  F  •  are  given,  following  the  same  procedure  as 

nj 

for  E^,  by: 

Fnj  *  WniO(-l)n[ft/(fa  *  ft)l  (C'30) 

where  f.  and  f*  are  given  by  Equations  C-59  and  C-60,  in  which  m  is 

Q  w 

2  2 

replaced  by  n(  *  1,  2,  3,  ...)  and  p  ■  ,  instead  of  P  *  (1T/ 2 } 

shown  after  Equation  C-60. 


114 


CASE  (3).  Constant  Heat  Input  at  One  End  and 
Insulated  at  the  Other 

The  problem  of  a  constant  heat  input  at  one  surface  while  the 
other  is  perfectly  insulated,  is  of  fundamental  importance  in  assessing 
different  rates  of  the  surface  temperature  rises  in  the  two  different 
regions.  Additionaly,  the  rates  of  heat  propagation  are  inter-mingled 
with  transverse  conduction  between  the  two  regions.  To  analyze  the 
temperature  response,  an  equivalent  thermal  capacity  (pc)e  is  defined 
on  a  prorata  basis  as: 

(pc)e  3  [(pc^a  +  (pc)2(b  -  a)]/b  (C-31) 

which  can  be  expressed  in  terms  of  the  non-dimensional  quantities  as: 

(pc)e  3  [(pc),.  +  (pc)2(b  -  a)]/b  (C-31 a ) 

Partial  Solutions.  Without  regard  to  the  interfacial  heat 
conduction,  a  "steady-state"  temperature  response  for  both  regions 
can  be  envisioned:  the  temperatures  will  be  linearly  (in  time) 
increasing.  Writing  the  partial  solutions  as: 


115 


(C-32a) 


T1  *  tQ/(2klx)H2/L3  +  [Q6/{pc)eL]  +  7, 

T2  *  [Q/(2k2x)][2/L]  +  [Q0/(pc)eL]  +  Yz  (C-32b) 

The  first  terms  of  the  preceding  equations  are  to  comply  with  the 
constant  heat  flux  Q  at  X  3  L,  and  the  second  terms  indicate  that 
both  regions  show  a  linear  temperature  rise  (equal  in  both  regions) 

In  time  as  a  consequence  of  the  heat  input. 

The  functions  71  and  ?2  are  to  satisfy  their  respective  governing 
equations.  Thus,  by  direct  substitutions  into  Equations  C-12  and 
C-13,  the  forms  of  and  Y2  can  be  determined,  and  Equations  C-32a 
and  C-33b  can  be  re-cast  as: 


3  (1^/2)  +  [(pc)1/(pc)e]9  +  [(pcyfpcjg  -  l](y2/2)/X]  (C-33a) 

T2  3  (kix/k2x^/2)  +  [(^/(pc^ie 

+(klx/k2x)t(pc)2/(pc)e  “  UC(y  “  b)2/2J/X2  (C-33b) 

where  Tj  3  ^kixTl^L^  and  T2  ■  (k]XT2/QL)*  (C-34) 

To  each  of  Equations  C-33a  and  C-33b,  a  constant  term  and  a  general 
solution  of  the  Laplace  equation  can  be  added.  The  potential  solutions 
from  their  respective  Laplace  equations  are: 
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cos(mrx)  cosh(mry/A^) 


and  cos(mrx)  cosh[mr(y  -  respectively. 


The  expanded  versions  of  the  partial  solutions  then  assume  the  forms: 


r  *  i/72  iu  1 
T1  I(x  '  3}  +  2x7 


7^  -  ']  ^  '  r  2)  *  f(«>,/(Pc>eJe 


+  £  Gi„  cos  ( mrx )  [cosh ( rnry/  A7)/cosh ( mra/  AT) ]  (C-35a) 

n*1 ,2,3,. .  ,n  1  1 


T2  *  (klx/k2x) 


1/72  h  +  1 

i(x  *  v  +  2xJ 


Rpc)?  1  _ 


E)2  -  (a  -  E)2] 


_  •  cosh[mr(y  -  F)/AT] 

+  [(pc)|/(oc)  10  +  £  c?„  cos(nirx)  - — -  (C-35b) 

nal,2,3,..  cosh[ror(a  -  F)/AT] 


The  coefficients  Gln  and  G2n  are  connected  by  the  condition  of  equal  heat 
flux  at  the  interface  y  s  a  through  the  relation: 


G1b  *  (k9u/ki„)ttan(mr(a  -  b)/ AT)/tan(mFa/.fi7)]G.,n 


(C-36) 


Next,  the  requirement  of  equal  temperature  at  y  ■  ?  is  imposed. 

Setting  y  *  T  in  Equations  C-35a  and  C-35b,  and  Tj  *  T2,  there 
results: 

(1/2)11  -(k,,/k?J)l^  -  1/31  *  £  ,  (5,„  -  GjJcosfmtx)  *  0 

(C-37) 

Through  a  regular  Fourier  series  technique,  the  coefficients  are 
obtained  as  follows: 


J2n 


2(1-  (k1v/k2<)l(-l)n 


2  2 
n  w 


1  -  <kiy/k2y) 


tan(mt(?  -  b)/A^) 
tan(niraT/A^) 


(C-38) 


The  partial  solutions  representing  the  steady-state  behaviors  are 
now  defined  by  Equations  C-35a,  C-35b,  C-36  and  C-38 


Transient  Part.  Since  the  Initial  conditions  of  *  T2  *  0 
are  yet  unfulfilled,  transient  solutions  of  the  following  forms  are 
needed : 


T,  »  £  E  H  .e  nJ  cos(mrx)^i  ni 

1  n*0, 1,2,..  j-1,2,..  nJ  l,nJ 


r2-  £  £  H  ,e  nJ  cos(mrx)72  . 

2  n-0,1,2,..  j-1,2,.  nJ 


(C-39a) 


(C-39b) 
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where. 


cosyjie^  -  (nir)2]A1 
cosa^/[8^  -  (mr)2]/x1 


(C-40a) 


and 


cos(y  -  b)yj [(o^/a^e^  -  (mr)2]A2 
cos  (a  -  Fj^Ka,/ a2)sjj  -  (nir)2]/X2 


(C-40b) 


Again,  the  forms  of  Equations  C-39  and  C-40  assure  equality  of 
temperatures  Tj  and  To  at  y  *  a\  The  requirement  of  equal  heat 
flux  results  in  the  following  eigen-value  equation  for  6nj ( j  *  1,2,..) 
for  each  n-value: 


-  (nir)2]/X1  tana  J[0^  -  (mr)2^ 


+  (k2y/kly^^°i/a2^Bnj  "  (nrr)2 J/Xg  tan(b  -  a)  yj  [(cyagjBjjj  -  (nn)2]/X2  *  0 


(C-4T ) 


for  n  *  0,1,2,  ...  Note  that  Equation  C-41  is  identical  to 
Equation  C-29. 
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The  complete  solutions  for  these  two  regions  are  therefore: 

T|  ■  Equation  C-35a  +  Equation  C-39a  (C-42a) 

T2  •  Equation  C-35b  +  Equation  C-39b  (C-42b) 

together  with  Equations  C-36,  C-38,  C-40a,  C-40b  and  Equation  C-41 

as  auxiliary  equations  defining  the  various  parameters.  The  coefficients 

H  .  are  evaluated  in  Section  C-III. 
nj 


CASE  (4).  Constant  Heat  Flux  at  One  End  and  Initial 
Temperature  Maintained  at  the  Other  End 


The  temperature  distributions  In  the  two  regions  are  expressed  by: 


oo 

Tj  =  x  +  £  Amsin(imrx/2)  [cosh(imry/2»^j')/cosh(nnra/2A^)] 

m»l ,3 )5| « • 


00  00  m  0 

+  £  L  s  .e  s1n(um7/2)7,  .  (C-43a) 

m-1,3,5,..  J-1.2...  ** 
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T2  «  (klx/k2x)x 


+  £  Bms1n(nrnx/2)tcosh(imr(y  -  F)/2/fp/cosh(imr(a  -  F)/2/fp] 

hi*  1  i3i5t*  • 


+  5Z  £  S  ,e  s1n(rmrx/2)7,  .  (C-43b) 

m*l,3,5,. .  j*l ,2, . .  Z_mJ 

The  eigen-functions  ^  and  ?2  ^  are  as  follows: 

71-mj  *  -  (imr/2)2]/X1  j cosa  -  (imr/2)2]/A1  j  (C-44a) 

Y2-mj  *  Top/Bottom  (C-44b) 

where  Top  •  cos (y  -  F) 

l(ol/a2)emj  '  ^2)2]/A2 

and  Bottom  ■  cos(?  -  F)  ^[(a.|  /a2 ) 6mj  ’  (mw/2)2]/X2 

The  eigen-values  are  defined  by  the  roots  of  Equation  C-21. 

In  equations  C-43a  and  C-43b,  the  coefficients  ^  and  Bffl  are  obtained 
by  temperature  and  heat  flux  equalities  at  y  ■  a.  Heat  flux  equality 
yields: 
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Am  *  "Bm(k2yy,,c1y^/Al/X2^tanh^'mr^  '  a)/2Ap/tanh(mira/2Ap]  (C-45) 

By  setting  T^(y  ■  ?)  ■  TgCy  ■  T),  there  ensues  from  Equations  C-43a 
and  C-43b  (for  x  ■  0  to  1), 


x[l  -  (klx/k2x)]  -  £  <Bm  -  Am)s1n(imrx/2)  (C-46) 


Since  x  has  a  Fourier  expansion  of, 

.  m-1 

x  *  L  l8/(imr)2](-l)'7’s1n(m1r7/2) 

m® 1 ,3,5, « i 


hence. 


B 


m 


(C-47) 


Equations  C-45  and  C-47  therefore  are  used  to  obtain  Am  and  Bm  as 
follows: 


r  m-n 


A.  tanh(mrr(b  -  a)/2A^) 

—  i  ,n  — — ■■  i 

tanh(imra/2VXp 


(C-48) 


122 


Equations  C-45  and  C-48  therefore  fix  the  coefficients  Am  and  Bm<  The 
double  Fourier  coefficients  Smj  in  Equations  C-43a  and  C-43b  are 
derived  in  Section  C-II - 


C-II.  SEGMENTAL  ORTHOGONALITY  OF  Y. 

mj 


The  function  Yj^  is  defined  by  the  following  expressions  and  is 
governed  by  two  differential  equations  in  two  sub-regions. 


Tmj  2,mj 


T  _  C0$y> 

/  <4l  ■ 

l»mj 

cosa^ 

Kj  -  p"2IAl 

a*,  and 

r 

cos(y  -  b)yj 

«YYsri  -  p">2>a2 

cos (a  -  F) 

> 

[(^/Og)^  -  pm2]/x2 

(C-49) 


(C-50) 


from  y  ■  T  to  F.  In  both  expressions,  p  is  simply  a  parameter. 

2 

Thus,  to  bring  about  Equations  C-20a  and  C-20b,  p  is  set  to  it  /4. 
In  what  follows,  m  does  not  need  to  be  an  integer. 

Equations  C-49  and  C-50  are  solutions  to  their  respective 
differential  equations  given  on  the  next  page. 
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xl7i,l«j  *  I<r  p"2|71,«J  *  0  (for  y  -  0  to  a)  (C-51) 

x272.mj  *[(V«2lBmJ  *  ’  °  (for  y  -  ?  to  F)  (C-52) 


Multiply  Equation  C-51  by  7^  mi  and  Integrate  both  sides  from  y  «  0 
to  a.  Similarly,  multiply  Equation  C-52  by  (k2X/klx^2  mi  and 
Integrate  both  sides  from  y  a  7  to  F.  The  results  are  (re-calling 

71,.U(I)  *72.«j®  ■ 

a  7 

xi7i.*®  '  xi  /  Ti  ,mj7i  '  i<j-  P"2'/7,.^7!,-^  <c-53> 


0 

-A2(k2x/klx)Y^mj(7)  -  X2^k2x/klx^/72,mj72,m1dy 


«  |(«1/c»2)(k2x/klx)0^j  -  (k2x/klxJpinZ]JL  Y2,mjY2,m1dy  (C'54) 

7 

In  obtaining  the  above  expressions,  the  fact  that  Yj  ^ (0)  *  0,  and 
Y2  mj^  *  0  has  been  taken  Into  account.  Furthermore,  It  can  be 
shown  that: 


Xl7l.mj^  "X2(k2x/klx^2,mj(a) 
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Is  identical  to  the  eigen-equation  of  heat  flux  equality  at  the 
interface,  y  *  a.  By  adding  Equations  C-53  and  C-54,  there  results: 


t?  b 

o  I 


a  o 

Xl/Yl.m/i,midy  “  X2(k2x/kix)/ 72,mjY2,midy 


(C-55) 


By  interchanging  the  Indexes  1  and  j,  there  ensues  another  equation 

which  has  the  right-hand-side  identical  to  that  of  Equation  C-55.  On 

the  left-hand-side,  however,  the  only  difference  is  6^2  Instead  of 
2 

8,^  •  Since  these  two  eigen-values  are  distinct,  hence  it  follows: 


/7l,m7l,«1d>r  +  /«» cy<sc)i  l^illv)72iBtdy  - 


(C-56) 


Thus,  the  segmental  orthogonality  of  Y^'s  is  established  with  respect 
to  the  segmental  weight  function  defined  by  Equation  C-22. 


C.III.  THE  COEFFICIENTS  E^  IN  EQUATIONS  C-U  AND  C-15 


Combining  Equations  C-14  and  C-15  Into  a  single  equation, 
these  results: 


T 


QO 

1  ♦  z 


-0  20" 

Z  Em<e  cos(imrx/2)7  . 

m\ ,3,5,. .  ^ 


(C-57) 


where  are  those  defined  by  Equations  C-20a  and  C-20b  in  their 
respective  regions.  In  order  to  satisfy  the  initial  condition  of 
T  *  0  at  9  *  0,  Equation  C-57  becomes: 


00  00 


j«l,2,..  m*l,3,5. 


E«jcos(”r/2)?«j 


-1 


multiplying  both  sides  by  cos(imrx/2)Ym^wf  and  Integrate  the  resulting 
products  on  both  sides  from  3T  ■  0  to  1  and  y  ■  0  to  b,  the  result  Is: 


E 


nl 


-4 

TmFT 


In  the  procedure  above,  the  orthogonality  relations  of: 


1 


cos  (imix/2)cos  (nnx/2  )dx 


0 
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I 


for  m  f  n,  and 


U 


for  i  f  j  have  been  incorporated.  The  complete  expression  for  E  .  is 

mj 

therefore: 


Emj  *  -(8/"w)(-l)  2  fft/(fa  +  ft)] 


(C-58) 


ft  -  I  tana  J[(*J  -  (pm2)]/X1  |  /J[ S^2  -  (pm2)]/*, 


+  [(oc)2/(pc)1l  tan(¥  -  a)^  Koj/og)^2  -  ^l/XgJ^K^/og)^2  -  pm2]/X2 


(C-59) 


f a  *  a  /  Jcosa  -  pm2]/^ 

+  [(pc)2/(pc)1 1 C5”  -  a)  j | cos(F  -  D^lfaj/agJe^,2  -  pm2]A2  ... 


(C-60) 


2 

where  p  *  (ir/2)  . 
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C.IV.  ANALYSIS  OF  EIGEN-VALUE  Bm,  AND  8  . 

mj  nj 


The  eigen-value  Equation  C-21  applies  to  Case  1  and  Case  4. 
and  the  eigen-value  Equation  C-29  applies  to  Case  2  and  Case  3.  These 
two  Equations  C-21  and  C-29  can  be  lumped  together  as: 


<Sf> 


] /X i  tana 


-  (if)‘ 


1/x, 


+  (k2y/kly)/i(al/ot2)BijZ  "  (t)_1A2  tan(¥  '  «)  J«Va2)6ij2  ’  = 


tln.y 


(C-61) 


For  1  «  1,3,5,  ....  then  the  above  equation  reverts  to  Equation  C-21; 
for  1  ■  0,2, 4, 6,  ....  the  above  equation  reverts  to  Equation  C-29. 

For  all  four  cases,  the  eigen-values  B.jj(j  3  1,2,  ...)  are,  of  course, 
all  positive  and  real. 

In  order  to  facilitate  a  unified  numerical  procedure  in  obtaining 
the  roots  B^j,  two  separate  cases  are  considered:  (i)  (a.,/a2)  >  1 
and  (11)  (a-i/c^,)  <  1. 

(1)  (c^/ag)  >1:  Let  4»a  ■  (c^/oig). 
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together  with  the  following  substitutions: 

*k  “  <k2y/k1y) 

^  =  1 

<|>2  *  (5/1)  -  1 

(V>  =  (2/ir)V 

9r  xi 

g2  “X2 


then  Equation  C-61  becomes: 

V KSija)2  -  (ia)2l/gi  tan  fc,  ^[(I^a)2  -  (ia)2 

+  <l>k\ll<t>aC& i/)2  -  Ua)2]/g2  tan(^2yj  t^C^ja)2 


(i  =  0,  1,  2,  3,  ...) 


(11)  (a-j/a2)  <  1: 


Let  <|>a  =  (o2/a^) 


together  the  following  substitutions: 

0k  "  (kly/k2y) 

<(»1  »  (F/a)  -  1 
$2  a  1 

(fT^a)  »  (2/tt) (0. jT)^(a1/a2) 

9l  s  x2 

g2»Xr 

then  Equation  C-61  becomes: 

\  >/[* aCSij-a)2  -  (ia)2]/g2  tan 

+  yj  [(Bjja)2  -  (15')2]/g1  (tan  tyiyji(  S^a)2  -  (ia)2]/g]  )  *  0  (C-63) 


Of  course.  Equations  C-62  and  C-63  are  identical,  which  Is  what  the 
two  sets  of  substitutions  are  designed  for. 

To  solve  for  (TT^F)  from  Equation  C-63,  first  consider  the  case 
for  1  f  0.  Since  $a  >  1,  hence  there  could  be  roots  of  (F^a)  such 
that  the  arguments  of  the  tangent  terms  in  Equation  C-63  become  of 
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opposite  signs.  The  roots  of  fl^a)  then  lie  between  the  lower  limit 
of  (Ta/,^]  and  the  upper  limit  of  (ta).  The  elgen-value  equation, 
reduced  from  Equation  C-63  becomes: 

-  yjl( la)2  -  ($ija)2]/g1  tanh^1  >/[(1a)2  -  (s^a)2]/^^ 

+  '  (ia)2l/g2  tan  ^g^t^d^a)2  -  (ia)Z]/g2j  =  0 

(C-64) 

In  Equation  C-64,  the  first  term  is  always  negative,  therefore 
the  value  of  (^j^)  satisfying  Equation  C-64  must  be  such  that  the 
angular  argument  In  the  second  term  must  be  In  odd  quadrants.  Hence, 
a  straightforward  search  for  the  roots  can  be  initiated.  Of  course, 
when  1  a  0,  Equation  C-64  1$  moot  and  the  full  Equation  C-62  must  be 
used. 

For  the  succeeding  range,  l.e.,  fiT^a-)  >(1T),  solution  for  the 
roots  can  be  simplified  by  observing  that  the  terms  in  Equation  C-62 
must  be  of  opposite  signs.  Hence,  by  starting  with  the  argument  of  the 
first  term  In  the  first  quadrant,  then  the  equation  is  satisfied  by 
the  argument  of  the  second  term  located  In  the  second  or  fourth 
quadrant.  The  procedure  Is  systematized  by  using  a  search  Index  K  for 
each  range  of  search  for  the  root. 


(i)  K  ■  0.  The  eigen-values  (fT^a)  are  determined  from  Equation  C-64 
and  are  bounded  by  the  following  limits: 


(1a)2/<J>a  <  (I^a)2  <  (la)2 


The  argument  in  the  second  term  of  Equation  C-64  must  be  in  an 

odd  quadrant  to  satisfy  the  equation. 

Of  course,  the  search  described  In  this  range  is  only 
valid  if  i  t  0.  Otherwise,  the  search  must  start  with  K  =  1 . 


(ii)  K  *  1.  Here  the  index  K  a  1  indicates  that  search  for  the  roots 
(B^a)  such  that  the  argument  of  the  first  term  of  Equation  C-62 
Is  In  the  first  quadrant.  Accordingly,  the  argument  of  the 
second  term  of  Equation  C-62  must  be  in  an  even  quadrant.  The 
range  of  (B^a)  can  be  expressed  by: 


(K  -  l)f  <  V((*tja)2  -  (1a)2]/9l  <  Kj 


(g1  (K  -  1  )2/*f  +  0a)2]  <  (S^a)2  <  (K/^ )2  +  (la)2] 


Within  these  two  limits  of  (TT^a),  possible  roots  of  (IT^a) 
are  to  be  located  such  that  the  argument  of  the  second  term  is 
in  the  second  or  fourth  quadrants. 
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(Ill)  K  »  2.  In  this  range,  the  argument  of  the  first  term  is  in 
the  second  quadrant  and  the  lower  and  upper  limits  on 
(B^T)  are  defined  by  those  shown  for  K  *  1.  The  roots  of 
O^a-)  are  those  for  which  the  argument  of  the  second  term 
Equation  C-62  is  in  an  odd  quadrant. 
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